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I.  INTRODUCTION 


It  is  frequently  important  for  the  analysis  of  propagation  of 
sound  in  the  ocean  to  be  able  to  compute  the  acoustic  field  in  a  range 
dependent  ocean  environment.  Although  good  computational  methods  have 
been  developed  for  modeling  acoustic  propagation  in  horizontally 
stratified  acoustic  media,  range  variations  in  the  environment  play  such 
an  Important  role  in  many  realistic  propagation  problems  that  the 
simplifying  assumption  of  horizontal  stratification  is  not  justified. 

The  methods  available  for  modeling  sound  propagation  in  a  range 
variable  medium  may  be  roughly  categorized  as  wave  theoretical  and  ray 
theoretical.  Reference  1  is  a  survey  of  these  methods. 

Wave  theoretical  methods  applicable  to  range  dependent 
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environments  include  the  adiabatic  normal  mode  approach,  horizontal  ray 
theory,3*4  the  parabolic  equation  method,5  coupled  mode  theory,6,7  and 
direct  numerical  solution  of  the  wave  equation.  All  but  the  last  of 
these  procedures  make  use  of  various  approximations  in  order  to  reduce 
the  required  computational  effort.  The  approximations  are  valid  only  for 
gradual  variations  of  the  media  with  range.  The  direct  numerical 
approach  is  computationally  very  intensive  and  may  not  be  practical  to 
Implement  except  on  the  largest  and  fastest  scientific  computers.  For 
all  of  these  models,  the  amount  of  computation  required  increases  rapidly 
with  Increasing  acoustic  wave  number,  so  the  models  are  effectively 
limited  to  low  frequency  applications.  The  principal  advantages  of  wave 
models  are  the  correct  treatment  of  diffraction  effects  and  the  ability 
to  correctly  model  complex  bottom  interaction  effects. 

The  principal  advantage  of  ray  theory,  on  the  other  hand,  is 
that  the  paths  of  energy  propagation  through  the  medium  are  explicitly 
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identified.  The  ability  to  single  out  the  energetically  important 
propagation  paths  appeals  to  physical  intuition  and  is  very  useful  for 
gaining  an  understanding  of  the  important  influences  on  propagation. 

The  convenience  of  the  ray  theory  formulation  is  achieved  at 
the  cost  of  an  approximation.  The  ray  theory  approximation  and  the 
approximations  made  in  practical  wave  models,  together  with  the 
computational  limitations  on  wave  models,  are  such  that  the  advantages 
and  limitations  of  ray  models  and  wave  models  are  almost  the  reverse  of 
each  other.  Wave  models  are  limited  to  low  frequencies;  ray  models  are 
most  accurate  at  high  frequencies.  Wave  models  include  diffraction  and 
boundary  interaction  effects;  uncorrected  ray  models  neglect  diffraction 
effects  and  have  less  sophisticated  treatments  of  boundaries.  Wave 
models  are  inaccurate  to  varying  degrees  near  sources,  particularly  in 
the  extreme  nearfield;  ray  models  grow  more  accurate  as  one  approaches 
the  source.  Most  wave  models  can  tolerate  almost  any  variation  of  the 
environment  with  depth,  but  ray  models  can  better  tolerate  range 
variations. 

In  recent  years,  ray  models  have  benefited  from  the  development 
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of  correction  formulas  which  introduce  diffraction  effects  ’  and  improve 
the  treatment  of  boundaries.*0 

The  author's  involvement  with  ray  models  began  in  1976  with  his 

development  and  implementation  of  a  range  independent  ray  model1*  at 

Applied  Research  Laboratories,  The  University  of  Texas  at  Austin 

(ARL:UT).  Prior  to  that  time,  several  range  dependent  ray  models  had 
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already  been  implemented  at  other  facilities.  ’  All  of  these,  however, 
suffered  from  various  deficiencies,  including  the  representation  of  the 
acoustic  media,  methods  for  tracing  rays  (determining  the  paths  of  energy 
propagation),  the  location  of  eigenrays  (finding  the  rays  which  travel 
from  the  source  to  an  observation  point  or  receiver),  the  calculation  of 
acoustic  pressure  along  a  ray  path,  and  the  design  and  implementation  of 
the  software  Itself.  These  problems  are  not  at  all  independent  of  each 
other  and  the  adverse  consequences  of  an  inadequate  solution  to  any  one 


of  them  will  usually  be  compounded  in  several  others.  Nor  could  these 
difficulties  reasonably  be  anticipated;  all  were  uncovered  by  attempting 
to  build  and  use  practical  models. 

In  1980,  with  U.S.  Navy  funding,  the  author  undertook  to 
develop  a  new  range  dependent  ray  model.  The  effort  drew  on  several 
years  of  experience  in  numerical  modeling  and  borrowed  heavily  from 
existing  ray  models  and  published  research.  The  result  is  the  ray  model 
MEDUSA,  which,  after  three  years  of  development,  testing,  and  revision, 
is  being  released  to  a  growing  number  of  research  and  prediction 
facilities. 

The  methods  developed  to  redress  the  problems  inherent  in  ray 
modeling  are  the  subject  of  this  report.  The  techniques  are  developed 
and  explained.  Considerable  attention  is  also  given  to  the  pitfalls  of 
ray  modeling  and  to  seemingly  plausible  procedures,  often  used,  which  do 
not  work.  Always,  the  focus  is  on  the  means  of  implementing  practical, 
reliable  ray  models  on  scientific  computers. 

The  remainder  of  this  report  is  organized  as  follows.  The  ray 
path  and  intensity  equations  are  derived  in  Chapter  II.  The  procedures 
for  solving  these  equations  are  described  in  Chapter  III.  The  methods 
for  computing  the  sound  speed,  sound  speed  gradients  with  respect  to 
depth  and  range,  and  the  bathymetry  (the  ocean  bottom  depth  as  a  function 
of  range  from  the  source),  given  only  tabulated  samples  of  these 
quantities,  are  described  in  Chapter  IV.  Chapter  V  covers  the  methods 
used  to  determine  the  eigenray  launch  angles  and  to  compute  eigenray  path 
data,  such  as  travel  times  and  arrival  angles.  Chapter  VI  is  an  overview 
of  the  organization  of  the  MEDUSA  model.  The  report  concludes  in 
Chapter  VII  with  a  typical  application  of  MEDUSA. 
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II.  REVIEW  OF  RAY  THEORY 


Brief  derivations  of  the  ray  path  and  intensity  equations  are 
presented  here  to  establish  the  notation  and  coordinate  system 
conventions  and  to  provide  background  for  the  development  of  the 
specialized  equations  used  in  MEDUSA.  More  extensive  treatments  of  ray 
theory  are  available  in  several  texts. 

The  acoustic  wave  equation  for  pressure  P(r,t), 


V2p  =  -Ty- - 

c  (r)  at2 


(II. 1) 


becomes,  upon  assuming  P(r,t)=A(r)e^^r^"^t, 


+  J 


2VA»V<(>  +  AV  <j> 


=  0 


(II. 2) 


where  k=u/c  is  the  wave  number  and  4>  represents  an  expanding  wave  front. 
Since  A  and  4>  are  assumed  to  be  real  functions,  the  terms  in  brackets  in 
Eq.  II. 2  must  vanish  separately  and  identically: 

l^l2  »  K2  ,  (II. 3) 

2$A-$<p  +  AV2<p  =  0  ,  (II.  4) 
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where  K^=k^+  V^A/A. 

Surfaces  of  constant  phase,  4)  ,  are  wave  fronts.  The  outward 
normals  to  these  surfaces,  parallel  to  the  vector  v<p  ,  indicate  the 
directions  of  wave  propagation.  Equation  II. 3,  the  eikonal  equation, 
gives  rise  to  equations  for  ray  paths  which  are  perpendicular  to  the  wave 
fronts.  Equation  II. 4,  the  transport  equation,  leads  to  equations  for 
intensity  along  the  ray  paths. 

The  intensity  at  a  point  on  a  ray  path  is16 


I  = 


* 

PP 


2pc  2pc 


(II. 5) 


where  p  is  the  density  of  the  medium.  For  two  points  located  on  a  ray 
and  both  points  in  the  water,  p  and  c  are  constant  to  within  a  few 
percent,  so  the  intensity  ratio  at  the  two  points  is  given  approximately 
by 


(II. 6) 


where  the  subscripts  indicate  the  two  points.  We  now  seek  an  expression 
for  A^/Ap 


The  volume  V,  shown  in  Fig.  1,  is  enclosed  by  the  surface  S, 
defined  by  the  wavefront  surface  sections  and  W£  and  the  four  ray 
paths  resulting  from  a  small  change  69$  in  the  ray  source  depression 
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angle  e$,  and  a  small  change  6^s  in  the  azimuthal  source  angle  ij>  . 
Noting  that  Eq.  II. 4  can  be  written  as  a  divergence, 


V*(A2^4>)  =  0 


(II. 7) 


one  obtains  by  Gauss'  theorem 


j ^  dV^-(A2^)  =J  dSn-A2^, 


(II. 8) 


where  n  is  the  outward  unit  vector  normal  to  S.  Since  n*V$=0  along  the 
ray  paths  by  definition,  the  only  contributions  to  the  integral  occur  at 
the  end  surfaces  and  W2.  Since  ^  is  always  directed  away  from  the 
source  one  has,  from  Eq.  I I. 3,  • 


A 

n 


(%)/K  on  Wj 
-(^4>)/K  on  W2 


(II. 9) 


Hence, 


(II. 10) 
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For  infinitesimal  variations  in  launch  angles,  the  integrands  will  be 
constants  on  the  surfaces  and  W2,  so  from  Eq.  II. 10, 


2  f 

a^k2  J w1 


(ii. U) 


Thus  far  the  only  approximation  made  has  been  the  assumption  that  p  and  c 
are  nearly  constant  in  the  water. 

If  A  varies  slowly  over  a  wavelength,  which  will  ordinarily  be 
true  for  high  frequencies,  slow  spatial  variations  in  c,  and  in  the 
absence  of  nearby  boundaries,  then 


k2  »  th 

K  »  A 


(11.12) 


Inequality  11.12  is  the  ray  theory  approximation.  When  it  is  valid,  K=k 
and,  since  k^k2  in  the  water. 


2  L  * 


h  i  •'“i 


(11.13) 


Equation  II. 13  states  that  the  intensity  ratio  is  inversely  proportional 
to  the  ratio  of  the  infinitesimal  areas  of  Wj  and  W 2* 
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We  now  seek  differential  equations  for  the  ray  paths.  Recall 
that  V<j  defines  the  direction  of  ray  propagation.  The  unit  vector  in  the 
direction  of  the  ray  path  is,  by  Eq.  11.3,  '.74>/K.  if  s  is  ray  path  length 
and  >T(s)  is  the  ray  path  position  vector,  then  dx/ds  is  also  the  unit 
vector  in  the  direction  of  the  ray  path,  so 


Differentiation  of  Eq.  11.14  with  respect  to  path  length  yields 


(11.14) 


(11.15) 


But,  using  the  vector  operator  identity  d/ds=dx/ds-7  and  the  fact  that 
V<t>=Kdx/ds  ,  the  right  hand  side  of  Eq.  11.15  can  be  simplified: 


<l  ?*  =  /*<.  fx  &  ♦ ;  at *  j  >t 

ds  ?  Ids  II  3x  y  3y  3z 


1  I  3*2 

2  K  7K 


(11.16) 
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If  we  again  invoke  the  ray  theory  approximation  (replacing  K  by  k),  then 
Eq.  11.16  becomes 


d_ 

ds 


$k 


(11.17) 


The  replacement  of  K  by  k  enables  one  to  factor  out  u  from  Eq.  11.17, 
making  explicit  the  frequency  independence  of  ray  theory.  Division  of 
Eq.  11.17  by  uj/c  yields 


(11.18) 


where  cQ  is  an  arbitrary  reference  sound  speed  and  n=cQ/c  is  the  index  of 
refraction.  Equation  11.18  is  the  ray  path  equation. 

The  two  basic  functions  of  a  ray  model  are  to  solve  the  ray 

path  equation  and  then  to  compute  intensities  at  desired  points  on  the 

ray  paths  by  determining  the  ratio  of  the  areas  of  and  W2>  In  the 
next  two  sections  differential  equations  are  derived  for  the  ray  paths 

and  intensities  in  forms  suitable  for  solution  on  a  digital  computer. 

A.  Ray  Path  Equations  in  an  Azimuthal ly  Symmetrical  Environment 

The  source  is  now  assumed  to  lie  on  the  vertical  axis  of  a 
cylindrical  coordinate  system  (see  Fig.  2)  and  the  sound  speed  and 
bathymetry  are  assumed  to  be  azimuthally  symmetrical  about  the  vertical 
axis. 
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1  1 


\ 


FIGURE  2 

RAY  PATH  SPREADING  IN  AN  AZIMUTHALLY  SYMMETRICAL  MEDIUM 


ARL:UT 
AS-82  1666 
TLF - GA 
10-19-82 
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If  we  equate  vector  components  on  the  left-  and  right-hand 
sides  of  the  ray  path  equation  (Eq.  11.18),  we  obtain  equations  for  the 
ray  path  position  coordinates  in  terms  of  the  path  length.  In  the 
cylindrical  coordinate  system  depicted  in  Fig.  2,  the  ray  path 
coordinates  are  range  (r),  depth  (z),  and  azimuthal  angle  (>//).  In  scalar 
form,  Eq.  11.18  becomes 


d_ 

ds 


in 

3z 


(II. 19a) 


and 


(II. 19b) 


(II. 20c) 


If,  for  convenience,  the  ray  azimuthal  angle  at  the  source,  4*  , 
is  taken  to  be  zero,  then  the  solution  to  Eq.  II. 19c  is  simply  ^ ( s )  =0, 
because  the  azimuthal  symmetry  of  the  environment  assures  that  rays  which 
originate  in  the  r-z  plane  remain  in  the  r-z  plane.  We  therefore  need 
only  concern  ourselves  with  Eqs.  II. 19a  and  II. 19b. 

Equations  II. 19a  and  II. 19b  are  parametric  second  order 
nonlinear  equations  for  the  ray  path  position  coordinates  r  and  z  in 
terms  of  the  path  length  s.  These  two  equations  can  be  combined  in  a 
single  second  order  nonlinear  differential  equation  for  path  depth  z  as  a 
function  of  range  r  by  changing  the  independent  path  parameter  from  path 
length  to  range.  The  change  of  variable  is  accomplished  using  the 
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relation  d/ds  =  cosed/dr,  where  0  is  the  vertical  ray  path  angle,  or 
depression  angle  (see  Fig.  2).  With  the  change  of  independent  variable, 
and  using  d^/ds=0,  Eqs.  1 1. 19a  and  1 1. 19b  become 


dz 

dr 


cosQ 


d_ 

dr 


(n  cose)  +  n  cos^e 


d2z 

d7 


n 


z 


and 


cose  (n  cose)  =  n,. 
dr  '  r 


where  nz  and  nr  are  the  depth  and  range  gradients  of  n. 
Eg.  II. 20b  into  Eq.  II. 20a  then  gives 


z" 


(l+z,2> 
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where  primes  denote  differentiation  with  respect  to  r. 

The  initial  values  required  for  solution  of  the 
are  obtained  from  the  source  position  and  ray  launch  angle: 


2(0)  - 


and 


z '  (0)  =  tane$ 


(II. 20a) 


(II. 20b) 


Substitution  of 


(11.21) 


path  equations 


(11.22) 


(11.23) 
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where  z  is  the  source  depth  and  0  is  the  launch  depression  angle  (see 
Fig.  2). 

Whenever  a  ray  undergoes  a  surface  or  bottom  reflection,  the 
solution  of  the  path  equation  must  be  restarted  at  the  boundary 
intersection  with  a  new  ray  path  slope.  The  ray  path  slope  upon  specular 
reflection,  z^,  is  calculated  from  the  incident  ray  path  slope,  zj,  and 
the  ocean  bottom  slope,  z^,  if  applicable.  In  the  case  of  a  surface 
reflection,  the  reflected  ray  path  slope  is 

=  -zf  .  (11.24) 

In  the  case  of  bottom  reflection. 


z£  *  tan(29g-9I) 


(11.25) 


where  9j=tan"^(zj)  is  the  incident  ray  path  angle  and  6g=tan"^(z^)  is  the 
ocean  bottom  slope  angle. 


B. 


Intensity  in  an  Azimuthally  Symmetrical  Environment 


As  shown  in  the  derivation  of  Eq.  11.13  and  in  Fig.  1,  the 
calculation  of  intensity  depends  upon  the  calculation  of  the 

infinitesimal  areas  of  the  wave  front  surfaces  Wj  and  W^,  which  are  swept 

out  by  making  small  variations  in  the  ray  launch  angles.  In  the 

azimuthally  symmetrical  environment  of  Fig.  2,  the  area  of  W2  at 

horizontal  range  r  is  given  by 


area  W2  =  (r6ip$)  (6z  cose)  = 


(r6if»s)  (fir  sine) 


Since  Eq.  11.13  is  an  expression  for  intensity  ratios,  it  is  convenient 
to  compute  the  ratio  of  the  intensity  at  an  arbitrary  point  along  the  ray 
path  to  a  reference  intensity  measured  one  unit  distance  from  the  source. 
The  area  of  the  surface  at  unit  distance  from  the  source,  W^,  is  from 
Eq.  11.26  and  Fig.  2, 


area 


COS9$  6ips  69$ 


(11.27) 


Then,  by  Eqs.  11.13,  11.26,  and  11.27,  the  intensity  ratio  is 
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r  cos9  6z/69$ 

r  sin9  6r/66s 

(11.28) 


Passing  to  the  limit  of  vanishing  launch  angle  variations. 
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(11.29) 


The  calculation  of  intensity  thus  requires  the  calculation  of  either 

3z/30$  lr-const  or  9r^"s 'z=const'  the  vertical  and  horizontal  rates  of 
ray  spreading  with  changing  launch  depression  angle. 

The  intensity  ratio  of  Eq.  11.29,  when  expressed  in  decibels, 
is  the  propagation  loss,  or  transmission  loss,  due  to  geometrical 
spreading  along  the  ray.  A  fuller  treatment  of  propagation  loss  must 
await  a  discussion  of  eigenrays,  and  so  is  deferred  to  Chapter  VII,  where 
the  role  and  Implementation  of  propagation  loss  in  ray  models  is 
described. 
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One  can  sometimes  estimate  3z/36$  by  tracing  two  rays  with 
slightly  different  launch  angles  out  to  range  r,  where  the  arrival  depths 
will  differ  by  6z  for  a  launch  angle  difference  of  60  .  One  then  makes 
the  finite  difference  approximation  9z/3©s -Sz/6©s .  Similarly,  one  might 
estimate  3r/39s  by  tracing  two  rays  to  a  common  depth  and  using 
3r/30s~6r/60s.  Indeed,  these  procedures  are  used  in  ray  models,  but  they 
are  not  very  satisfactory.  They  can  at  best  attain  an  error  on  the  order 
of  S0£  and  are  subject  to  cancellation  errors. 

Moreover,  even  a  slight  change  in  the  ray  launch  angle  may 
cause  a  ray  to  follow  a  path  which  suddenly  veers  from  the  course 
followed  by  the  companion  ray.  For  example,  a  slight  change  in  the 
launch  angle  of  a  ray  which  passes  near  the  surface  of  the  ocean  could 
cause  the  ray  tostrike_thg.surface_  .and-  ref4®efc  -fnom-4t—  --The— surface  - 
reflected  ray  and  the  surface  grazing  ray  would  not  be  suitable  ray  pairs 
for  the  finite  difference  method  at  any  range  beyond  the  surface 
reflection  point  because  both  3z/3©s  and  3r/39$  are  discontinuous  in  the 
launch  angle  interval  (0  0  +60  )  beyond  that  range. 

To  avoid  numerical  differentiation,  we  seek  a  method  of 

computing  3z/30.  along  a  ray  path  as  the  ray  trace  progresses.  Solomon 
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and  Armijo  proposed  a  suitable  method,  described  below. 

The  ray  path  depth  z  at  range  r  is  a  function  not  only  of  the 
range,  but  also  of  the  ray  launch  angle,  that  is,  z=z(r,0$).  One  may 
therefore  differentiate  the  ray  path  equation  (Eq.  11.21)  with  respect  to 
0$  to  obtain 


C" 


2z'z 
1  +  z 
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(l+z'2) 


,  (11.30) 
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where  s=3z/39s,  nzz=3  n/3z  ,  nrz=3  r/3z9r,  and  primes  denote  d/dr. 

Initial  values  £  (z$)  and  ^'(z$)  are  needed  to  solve  Eq.  11.30. 
Since  the  source  depth  is  given  and  does  not  depend  on  0  , 


SUS)  =  0 


(11.31) 


Since  z' (r)=tane. 


9z 1 
36- 


36, 


2 

sec  0 


so  that,  at  the  source. 


?’(zs)  =  sec2es 


(11.32) 


(11.33) 


When  a  ray  reflects  from  the  bottom,  a  new  set  of  initial 
values  for  c  and  V  is  needed  in  order  to  proceed  with  the  solution  of 
Eq.  11.30.  Figure  3  shows  two  rays,  differing  in  launch  angle  by  the 
small  increment  fi6s,  striking  the  bottom  at  two  slightly  separated 
points,  labeled  1  and  2.  The  value  of  r,  at  the  point  of  incidence  is,  by 
definition,  the  limit 
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1  im 
66s-0 


(11.34) 


(see  Fig.  3). 
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Similarly,  the  value  of  4’  at  the  point  of  incidence  is 
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FIGURE  3 

THE  EFFECTS  OF  CHANGING  LAUNCH  ANGLE 
ON  BOTTOM  REFLECTED  RAYS 


(11.35) 


6e$-0  6e$ 


Upon  reflection,  5  and  r,'  take  on  the  new  values 


1  im 
66s->0 


and 


X 


R 


1  im 


(11.36) 


(11.37) 


where  z^  is  the  ray  path  slope  at  point  1  upon  reflection.  We  will  now 
derive  expressions  for  ?R  and  in  terms  of  Cj,  Cj,  and  the  geometry  of 
the  ray  path  and  bottom.  The  task  is  complicated  by  the  fact  that  both 
the  ray  paths  and  the  bathymetry  are  curvilinear. 

The  approach  used  will  be  to  express  the  values  of  z,  z',  X  , 
and  4'  at  the  points  2,  3,  and  4  as  first  order  Taylor  series  expansions 
about  point  1  and  substitute  in  Eqs.  11.36  and  11.37.  The  final 
expressions  are  exact,  not  first  order  approximations,  because  of  the 
limiting  process. 

Let  us  begin  with  the  derivation  of  by  deriving  an 

expression  for  z^  to  first  order  in  6e$. 
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(11.38) 


Z3  =  z2  *  (rrr2)2R2  + 
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(11.39) 


where  5rBar2_ri>  zbi  is  the  bottom  slope  at  point  1,  z,}2  is  the  ray  path 
slope  upon  reflection  from  point  2,  and  3rB/3e$  is  the  rate  at  which  the 
range  of  the  bottom  intercept  changes  with  launch  angle.  Note  that  in 
Eq.  11.39,  it  is  permissible  to  make  the  approximation  because 

only  first  order  terms  in  66$  are  being  retained.  Later,  during  the 
derivation  of  the  full  first  order  expansion  of  z^2  will  be  needed. 

The  quantity  8^g/9f)s  will  be  needed  several  times  in  the 
remainder  of  this  section  and  also  in  Chapter  V,  so  it  is  worthwhile  to 
develop  an  expression  for  it  here.  Again,  to  first  order, 
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But  another  expansion  for  z ^  is  also  available: 


Z4  =  Z1  +  <Vi  +  *• 


(11.41) 


Z2  =  Z4  +  66s?4  + 


z1  +  6rgz'  +  60S(C1  + 


•  )  + 


Z1  *  «Vi  *  * 


(II. 42) 


Combining  Eqs.  11.40  and  11.42  and  taking  the  limit  69 -*-0  yields 
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(11.43) 


Notice  the  use  of  the  first  order  approximation  in  Eq.  11.42. 

With  z^  and  9 r g/ 9©s  now  in  hand,  it  is  possible  to  evaluate  the 
right-hand  side  of  Eq.  11.36  to  obtain 


(11.44) 


IB _ ‘R 

Cr  '  ~Cl  zF^T 

If  desired,  Eq.  11.44  can  be  rewritten  with  the  aid  of  trigonometric 
identities  in  the  form 


cose. 

CR  =  "?I  cose^  *  (11.45) 


Turning  now  to  we  see  from  Eq.  11.37  that  a  first  order 
expansion  of  is  needed. 


z3  =  ZR2  *  *rrr2*ZR2  *  — 
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Here,  we  are  able  to  make  the  approximation 
must  develop  the  first  order  expansion  of  z 
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z2  =  z4  +  60s«4  +  *'• 


Using 
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(11.49) 


-R2  =  tan(29B2'92) 
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.  (11.50) 
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tan(x+e)  =  tan  x  +  e  sec  x  +  ...  , 


(11.51) 


tan_1(x+e)  =  tan-1x  +  e  - + 

1  +  xz 


(11.52) 
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and  retaining  only  terms  up  to  first  order  in  50$,  Eq.  11.50  becomes 
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(11.53) 
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We  may  now  evaluate  the  right-hand  side  of  Eq.  11.37,  using  Eqs.  11.43, 
11.46,  and  11.53,  to  obtain 
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New  values  for  ?  and  ?  '  are  also  required  when  a  ray  reflects 
from  the  surface.  Derivation  of  expressions  tor?,  and  proceeds  in 
the  same  way  as  for  bottom  reflections,  but  is  very  much  simplified  by 
the  assumption  that  the  surface,  unlike  the  bottom,  is  horizontal  and 
perfectly  flat.  The  results  are 


and 


sR  =  “Si 


(11.55) 


Sr  =  S!  z 


(11.56) 
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Surmian 


In  this  chapter  expressions  involving  differential  equations 
have  been  derived  for  ray  paths  and  intensities  in  azimuthal ly 
symmetrical  ocean  environments.  Also  derived  were  initial  values  needed 
to  start  the  solutions  of  the  differential  equations  at  the  source  and  to 
restart  the  solutions  when  rays  reflect  from  the  surface  and  bottom. 
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III.  SOLUTION  OF  THE  RAY  EQUATIONS 


I 


In  this  section  we  will  discuss  methods  for  solving  the  ray 
path  equation  (Eq.  11.21)  and  the  depth  spreading  equation  (Eq.  11.30), 
which  are  practical  for  implementation  on  a  digital  computer.  Since  the 
computational  effort  expended  tracing  rays  typically  exceeds  by  two 
orders  of  magnitude  the  efforts  required  in  any  other  single 
computational  task,  and  since  the  success  of  any  attempt  to  locate 

eigenrays  or  perform  any  other  further  processing  depends  critically  on 
high  quality  ray  traces,  the  ray  tracing  techniques  may  be  the  single 
most  important  facet  of  a  ray  model . 

Ray  trace  methods  based  on  analytical  solutions  of  the  ray 
equations  are  discussed  in  the  section  below.  Models  which  rely  entirely 
on  these  methods  have  generally  not  been  very  successful.  The 

difficulties  with  analytical  methods  are  set  forth  in  some  detail  not 
only  because  they  are  frequently  employed,  but  also  because  hybrid 

analytical/numerical  techniques  show  promise  of  substantially  improving 
future  ray  tracing  techniques.  Hybrid  methods  employ  numerical 

perturbation  procedures  to  refine  approximate  analytical  computations, 
but  the  analytical  expressions  used  must  be  extricated  from  the  pitfalls 
described  below. 

The  numerical  ray  tracing  methods  currently  employed  in  MEDUSA 
are  discussed  in  Section  III .B. 

A.  Analytic  Solution  of  the  Ray  Equations 

One  can  sometimes  solve  the  ray  path  equation,  the  depth 

spreading  equation,  and  even  the  equations  of  intersection  of  rays  with 
boundaries,  when  n(r,z)  and  the  bathymetry  are  particularly  simple 
functions.  For  example,  if  ^(1/n)  is  a  constant  vector,  then  the  ray 
paths  can  be  shown  to  be  circular  arcs. 
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Furthermore,  one  can  approximate  more  complicated  n's  to  any 
desired  degree  of  accuracy  by  breaking  up  the  r-z  plane  into  a  patchwork 
of  smaller  domains  and  fitting  n(r,z)  in  each  domain  with  simple 
(tractable)  functions.  The  ray  paths  are  then  constructed  by  solving  the 
ray  path  equation  analytically  within  each  domain  traversed  and 
connecting  the  path  segments.  One  could,  for  example,  break  up  the  r-z 
plane  into  triangles  within  which  the  ray  paths  are  circular  arcs. 

The  availability  of  analytic  solutions  to  the  ray  equations 
makes  this  procedure  seem  attractive;  and  it  is,  in  fact,  used  in  a 
number  of  ray  models.  However,  the  method  suffers  from  several 
drawbacks,  which  are  examined  below. 

Perhaps  the  most  serious  difficulty  has  to  do  with  seemingly 
innocuous  limitations  of  the  fitting  functions  used  to  approximate  n. 
When  n  is  approximated  with  some  of  the  simplest  and  most  convenient 
fitting  functions,  such  as  the  one  which  gives  circular  ray  arcs,  the 
resulting  gradients  nz  and  nr  are  generally  mismatched  at  the  boundaries 
of  the  small  domains.  These  discontinuities  in  turn  give  rise  to 
discontinuities  in  z(r,0s) |r=consf  (Recall  that  the  derivative  of  this 
function  with  respect  to  0 s  is  the  depth  spreading  function  c.)  In  other 
words,  if  the  fitting  functions  do  not  preserve  the  smoothness  of  n,  then 
the  ray  paths  do  not  change  smoothly  with  changing  launch  angle.  These 
erratic  ray  paths  produce  well  publicized  anomalies  in  propagation  loss 

IQ 

curves,  such  as  false  caustics  and  shadow  zones;  but  more  importantly, 
they  greatly  complicate  the  process  of  locating  the  eigenrays  (rays 
connecting  the  source  and  receiver)  which  are  needed  to  compute 
propagation  loss.  Location  of  eigenrays  and  the  ramifications  of 
discontinuities  are  discussed  in  Chapter  V. 

One  might  solve  this  problem  by  finding  fitting  functions  which 
are  flexible  enough  to  permit  the  gradients  of  n  to  be  matched  at  the 
domain  boundaries,  but  the  author  is  not  aware  of  any  such  functions 
which  do  not  introduce  nonphysical  contortions  in  n  and  yet  are  still 


simple  enough  in  form  to  allow  solution  of  the  ray  equations.  And  even 
if  such  functions  could  be  found,  it  turns  out  that  even  discontinuities 
in  the  second  derivatives  of  n  can  cause  anomalies  in  propagation  loss 
calculations. 

One  might  also  mitigate  the  problem  by  creating  smaller  and 
more  numerous  domains,  with  smaller  gradient  mismatches,  in  order  to 
better  approximate  n,  but  this  may  impose  unacceptable  demands  on 
computer  resources.  Yet,  experience  shows  that  ray  model  predictions  are 
often  sensitive  to  the  number  and  location  of  domains,  an  indication  that 
the  number  of  domains  typically  used  is  too  small  to  produce  stable, 
reliable  predictions. 

A  second  set  of  difficulties  relates  to  the  analytic  solutions 
to  the  ray  equations.  It  is  sometimes  assumed  that  solutions  of  the  path 
equations  in  the  form  of  combinations  of  algebraic  or  elementary 

transcendental  functions  automatically  guarantees  accurate  computation 
since  these  simple  functions  are  readily  available  on  computers  and  can 
be  evaluated  with  great  accuracy.  This  is  not  at  all  the  case;  the 
solutions  to  the  ray  equations  usually  will  have  to  be  recast  in 

mathematically  equivalent  forms  which  are  resistant  to  cancellation 
errors  and  other  numerical  difficulties.  To  do  this  requires  specialized 
knowledge  of  numerical  analysis  and  is  something  of  an  art  in  any  case. 
It  is  also  necessary  to  ensure  that  the  ray  formulas  do  not  break  down 
for  certain  limiting  cases.  The  author  knows  of  several  ray  models  which 
fail  when  n  is  a  constant. 

Lastly,  analytic  solutions  of  the  ray  equations  are  inherently 
dependent  on  the  form  of  the  function  n.  If  n  is  changed,  then  new 
solutions  to  the  path  equations  must  be  derived.  For  example,  if  a  ray 
model  were  developed  in  which  ^(1/n)  was  a  constant  vector,  then  the 
corresponding  ray  paths  would  be  circular  arcs.  If,  later  on,  one  wished 
to  change  n{r,z)  so  that  V(n  )  was  a  constant  vector  instead,  then  one 
would  also  have  to  revise  the  part  of  the  model  concerned  with  tracing 
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rays  so  that  the  circular  arcs  become  parabolic  arcs.  The  depth 
spreading  formulas,  boundary  intersection  formulas  (if  any),  travel  time 
formulas,  path  length  formulas,  and  any  other  ray  formulas  would  also 
have  to  be  altered.  The  necessary  modifications  to  the  model  clearly 
could  be  very  extensive,  yet  the  ability  to  change  n  readily  is  a  very 
valuable  asset  to  a  ray  model,  particularly  if  the  model  is  intended  to 
be  a  research  tool. 

B.  Numerical  Solution  of  the  Ray  Equations 

All  of  the  foregoing  difficulties  are  avoided  (or  at  least 
exchanged  for  a  different  set  of  problems)  if  the  ray  equations  are 
solved  numerically. 

Methods  for  solving  ordinary  differential  equations  have  been 
the  subject  of  vigorous  research  within  the  field  of  numerical  analysis 
for  several  decades  and  a  great  many  procedures  have  been  developed.  It 
will  therefore  not  be  necessary  to  develop  any  new  techniques,  but  rather 
to  identify  the  methods  best  suited  to  solving  the  ray  equations.  As  we 
shall  see,  the  methods  of  choice  turn  out  to  be  some  of  the  oldest  and 
simplest. 


The  most  outstanding  fact  about  ray  traces  from  the  standpoint 
of  computational  procedure  is  that  the  progress  of  a  ray  is  highly 
eventful.  It  may  undergo  surface  reflections  and  bottom  reflections,  it 
may  pass  through  caustics  (points  where  S  =0  and  the  ray  theory 
approximation  K-k  breaks  down),  it  can  encounter  receivers,  it  may  pass 
through  turning  points  (where  z'=0),  and  so  on.  It  is  essential  to 
diagnose  these  events  as  they  occur,  locate  them  precisely  and,  in  the 
case  of  boundary  reflections,  restart  the  ray  trace  with  new  initial 
conditions.  This  at  once  removes  from  consideration  various  global 
methods  for  solving  differential  equations  which  rely  on  variational 
principles  and  which  are  not  well  suited  to  discovering  and  responding  to 
local  disruptions. 
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We  shall  therefore  restrict  our  attention  to  incremental 
methods,  which  will  allow  us  to  advance  the  solution  to  the  ray  path 
equation  in  a  series  of  small  range  steps  leading  away  from  the  source. 
After  each  tentative  step,  an  exhaustive  battery  of  tests  will  be 
performed  to  locate  any  significant  events  which  may  have  occurred  within 
the  range  step.  In  response  to  some  events,  such  as  boundary  encounters, 
it  will  be  necessary  to  reject  the  step  and  adjust  the  range  step  size  so 
that  the  ray  path  segment  ends  exactly  where  the  event  occurs  in  order  to 
reinitialize  and  properly  resume  the  ray  trace  after  the  event. 

In  addition,  it  has  proven  highly  advantageous  to  be  able  to 
vary  the  step  size  in  response  to  changes  in  the  environment.  In  order 
to  maintain  desired  accuracy  in  the  ray  trace,  a  step  size  of  no  more 
than  a  few  meters  may  be  possible  as  a  ray  undergoes  rapid  undulations 
near  the  surface,  where  nz  may  vary  rapidly.  Yet,  step  sizes  of  a 
kilometer  or  more  may  be  permissible  for  the  same  ray  as  it  travels  in 
deep  water,  where  nz  is  nearly  constant. 

The  search  has  now  been  narrowed  to  methods  which  are  easily 

started  (and  restarted)  and  permit  the  step  size  to  vary  continuously. 

These  requirements  unfortunately  rule  out  the  multistep  and  variable 

order  methods  which  have  occupied  much  of  the  attention  of  numerical 

analysts  in  recent  years.  When  applicable,  these  methods  can  be 

extremely  accurate  and  efficient,  but  they  require  somewhat  complicated 

and  inconvenient  initialization  procedures,  they  can  only  double  or  halve 

their  step  sizes  conveniently,  and  they  lose  efficiency  when  interrupted 

?0 

after  every  step.  We  are  left  with  Runge-Kutta  methods  and  Taylor 
series  expansions. 

In  fact,  we  are  confined  to  Runge-Kutta  and  Taylor  series 
methods  of  low  order.  The  Taylor  series  expansion  for  the  ray  path  depth 
z  at  range  r+Ar,  given  the  path  depth  and  its  range  derivatives  at  range 
r,  is 
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z(r+Ar)  =  z(r)  +  Arz'(r)  +  ^  Ar2z"(r)  +  ^  Ar3z"'(r)  +  ... 


( HI  .1 ) 

The  accuracy  of  the  expansion  is  limited  by  the  number  of  continuous 
derivatives  of  z  there  are  in  the  interval  (r,r+Ar);  accuracy  increases 
as  more  terms  are  added  until  a  term  with  a  discontinuous  derivative  is 
encountered.  The  order  of  the  expansion  is  given  by  the  largest  power  of 
A r  included  in  the  series.  Runge-Kutta  methods  do  not  depend  explicitly 
on  the  derivatives  as  Taylor  series  do,  but  they  are  derived  from  Taylor 
series  expansions  and  are  subject  to  exactly  the  same  accuracy 
limitations.  The  order  of  a  given  Runge-Kutta  method  is  the  same  as  the 
order  of  the  Taylor  series  it  was  derived  from.  Repeated  differentiation 
of  the  ray  path  equation. 


z"  =  (1+z'2) 


(>■■■;) 


shows  that  the  higher  derivatives  of  z  remain  continuous  only  as  long  as 
the  gradients  of  n  remain  continuous.  A  practical  method  is  developed  in 
Chapter  IV  which  achieves  two  continuous  gradients  of  n  at  domain 
boundaries  (all  gradients  are  continuous  within  the  domains).  With  that 
limitation,  experience  indicates  that  greatest  efficiency  is  achieved 
with  methods  of  about  order  three.  Higher  order  methods  lose  efficiency 
crossing  domain  boundaries,  where  their  greater  computational  costs  are 
not  offset  by  increased  accuracy  and  the  ability  to  take  longer  range 
steps. 


Taylor  series  methods  are  effectively  limited  to  second  order 
because  of  the  inconvenience  of  computing  higher  derivatives  of  z  and  5 . 
Therefore,  the  method  currently  implemented  in  MEDUSA  is  a  third  order 
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Runge-Kutta  method,  specifically,  Kutta's  third  order  rule.  Using  this 


method,  the  ray  path  depth  and  ray  path  slope  at  r+Ar  are  computed  from 
z(r)  and  z'(r)  by  application  of  the  formulas 


z(r+Ar) 

z(r) 

+  — P 

z ' (r)  +  4zj  +  z^ 

6 

z' (r+Ar) 

z'(r) 

z " ( r )  +  4z^  +  z£ 

rr  +  Ar  1 

rl 

r  +  ~2 

Z1 

z ( r )  +  z ' (r) 

zi 

z'(r)  +  ^  z"  (r ,z(r) ,z‘  (r) j 

A 

V 

~r  +  Ar 

1 

Z2 

z( r)  +  Ar^2zj-z'  (r)j 

z2 

z' (r)  +  Ar 

2z'j  -  z "  (r , z ( r  1  z '  ( r )) j 

z2_ 

Z  (  *  Z2  9  ^2  )  J 

z"(r,z,z* ) 


(l+z'2> 


nz( r*z )  (  nr(r,z)  ^ 

n(r,zl  "  2  n(r,z)  / 


,  (III. 2) 


(III. 3) 


.  (III. 4) 


(III. 5) 


This  simple  Runge-Kutta  formula  lacks  a  procedure  for 
estimating  the  error  incurred  in  taking  a  step.  One  independent  means  of 


estimating  the  error  is  based  on  a  recasting  of  the  ray  path  equation  in 
the  form  of  an  integro-differential  equation: 

fr  i  2 

n*  cos^  =  n^  cos20q  +/  dr  .  (III. 6) 

Jr  0 

Equation  III. 6  is  a  generalized  form  of  Snell's  law  which  is  applicable 

when  n  is  a  function  of  range  as  well  as  depth.  In  a  horizontally 

stratified  (range  independent)  medium,  n  =0  and  Eq.  III. 6  reduces  to 
2  2  2  r 

Snell's  law:  n^cos61=  nQcos  0Q.  One  can  differentiate  Eq.  III. 6  with 

respect  to  r  and  recover  the  ray  path  equation.  If  one  lets  rg=r  and 

Tj=r+Ar,  takes  a  step  using  Eqs.  III. 2-5,  and  then  evaluates  the  left- 

and  right-hand  sides  of  Eq.  III. 6,  then  the  inevitable  discrepancy 

between  the  left-  and  right-hand  side  is  a  measure  of  the  numerical  error 

incurred  in  the  step.  This  procedure  does  not  yield  error  estimates  for 

the  ray  path  depth  and  slope;  it  is  a  measure  of  how  well  the  predicted 

2  2 

and  actual  values  of  n^os  0^  agree. 

The  error  incurred  in  one  step  using  Kutta's  third  order  rule 

4 

is  roughly  proportional  to  Ar  .  If  the  estimated  error  in  a  step  is 
unacceptably  large,  then  the  step  is  rejected  and  attempted  anew  using  a 
step  size  reduced  in  accordance  with  the  fourth  power  law.  On  the  other 
hand,  if  the  error  is  much  smaller  than  allowable,  one  can  increase  the 
step  size  for  the  next  step,  again  using  the  fourth  power  law  as  a  guide. 

The  efficacy  of  Kutta's  third  order  rule,  together  with  the 
error  estimator,  has  been  severely  tested  by  comparing  solutions  obtained 
using  these  methods  with  known  solutions  for  analytically  tractable  forms 
of  n.  Not  only  do  the  rays  arrive  within  a  few  millimeters  of  the 
correct  depths  across  range  intervals  of  thousands  of  kilometers,  but  the 
errors  themselves  are  usually  very  close  to  the  specified  tolerances. 
Kutta's  third  order  rule  is  also  used  to  solve  the  depth  spreading 
equation. 


The  error  estimation  procedure  calls  for  the  numerical 

evaluation  of  a  definite  integral.  The  integrity  of  the  error  estimation 

procedure  requires  a  numerical  integrator  of  higher  order  than  Kutta’s 

third  order  rule,  so  that  the  integrator  itself  does  not  become  an 

important  source  of  error.  A  suitable  integrator  is  the  Euler-Maclaurin 
22 

formula 


1  (x, “X- ) 

f ( x ) dx  =  \  0  [f(x0)+f(x1)j 


(xrx0)  1  (x.-x«)  (.X 

T2  [f,(x1)‘f,(xo)]  +  720  f  ^ 


(III. 7) 


XQ  <  V  <  Xj 

The  error  term  (x,-xn)^f (v)/720  compares  favorably  with  the  error  term 
of  the  more  familiar  Simpson's  rule  integrator.  The  Simpson's  rule 
error  term  is  (xj-xQ)5f ^ ( v)/90. 

The  Euler-Maclaurin  integrator  requires  not  only  the  values  of 
the  integrand  at  Xq  and  Xp  but  also  the  values  of  the  derivative  of  the 

integrand.  Thus,  inspection  of  Eq.  III. 6  reveals  that  the  range 

2 

derivative  of  8(n  )/3r  is  ^eded.  It  is 


d  Mn2)  =  (l_  +  3\  3(n2) 

cF  9r  \  9r  dzJ  9r 


4 


+  nn. 


rr 


z'(n  n  +nn  ) 
v  z  r  rz 


(III. 8) 
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The  author  has  also  experimented  with  Runge-Kutta-Fehlberg 
methods2^  for  solving  the  ray  equations.  They  proved  successful, 
offering  the  ability  to  estimate  errors  and  readily  adjust  the  step  size 
to  accommodate  error  tolerances  and  variations  in  the  environment,  but 
their  greater  complexity  and  a  tendency  to  balk  at  domain  interfaces  made 
them  less  advantageous  than  the  procedure  described  earlier.  Runge- 
Kutta-Fehlberg  methods  are  general  purpose  algorithms  and,  as  such, 
cannot  make  use  of  specialized  information  about  the  differential 
equations  they  may  be  required  to  solve,  as  we  have  done  here. 

It  is  expected  that  further  research  into  ray  tracing 
techniques  will  be  carried  out.  The  present  method  is  more  than 
sufficiently  accurate,  but  it  is  somewhat  slow. 

C.  Path  Length  and  Travel  Time  Integrals 

The  integrals  for  path  length  and  travel  time  are 


s  =  /  (l+z,2}1/2  dr  (III. 9) 

Jo 


and 


n ( 1+z ' 2 ) 1/2  dr 


(III. 10) 


The  most  popular  numerical  method  for  computing  such  integrals,  Simpson's 
rule,  is  inappropriate  here  because  the  range  points  at  which  the 
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1 


integrands  are  available  are  neither  equally  spaced  nor  necessarily  odd 
in  numoer,  as  required  by  Simpson's  rule.  However,  at  each  end  of  a 
range  step  one  has  available  the  integrands  and  their  range  derivatives 
at  both  ends  of  a  given  range  step: 


ds  _ 
dr 


(l+z'2)1/2  =  sec6 


(III. 11) 


^-5-  =  z’z"cos0 

dr 


(III. 12) 


dt  _  n  ds 

dF'^d?  *  (HI- 13) 

and 


(III. 14) 


The  third  order  Euler-Maclaurin  formula  is  used  to  compute  the 
contributions  to  the  integrals  over  a  range  step  using  the  integrands  and 
their  range  derivatives. 

D.  Summary 

In  this  chapter  the  ray  path  and  intensity  differential 
equations  were  cast  in  forms  suitable  for  solution  on  a  digital  computer. 
The  numerical  procedure  used,  Kutta's  third  order  rule,  was  described  and 
the  reasons  for  selecting  it  were  presented.  An  expression  for  the  ray 
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path  error  tolerance  was  derived  and  the  numerical  method  for  achieving 
the  accuracy  requirement  by  varying  the  integration  step  size  was 
presented.  Lastly,  the  use  of  the  Euler-Mac laurin  formula  to  perform  the 
ray  path  length  and  travel  time  integrations  was  explained. 
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IV.  REPRESENTATION  OF  THE  BATHYMETRY  AND  SOUND  SPEED 


Bathymetry  is  usually  presented  to  ray  models  in  the  form  of  a 
table  of  measured  ocean  bottom  depths  as  a  function  of  range  from  the 
source.  To  simplify  location  of  eigenrays,  it  is  necessary  to  fit  the 
tabular  data  with  a  function  which  provides  two  continuous  derivatives  of 
the  depth  with  respect  to  range.  There  is  a  function  which  is  designed 

pc 

to  do  just  that;  it  is  called  a  cubic  spline.  The  use  of  cubic  splines 
to  fit  sound  speed  profile  data  in  range  invariant  environments  was 

pc  1  O 

suggested  by  Moler  and  Solomon.  The  ray  model  GRASS  fitted  splines 
to  profiles  in  range  dependent  environments  and  traced  rays  numerically. 


Cubic  Splines 


The  cubic  spline  is  usually  defined  to  have  the  following 
properties: 


(1)  Given  tabular  data  x^,y^,  where  the  x^'s  are  abscissas, 
the  y^s  are  the  corresponding  ordinates,  and  l<i<n,  the 
cubic  spline  S(x)  is  defined  to  be  the  union  of  cubic 
polynomials  S^x).  Each  Si  is  defined  on  the  interval 
Xi<x<xi+r 

(2)  Si(xi)  =  yr 

(3)  Si (xi )  =  Sil(x.). 

(4)  S.(xi)  =  S._1(xi). 

(5)  SI<X1>  ■  ■  0  . 
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It  can  be  shown  that  the  cubic  spline,  thus  defined,  is  unique. 
It  can  also  be  shown  that 


(IV. 1) 


that  is,  the  cubic  spline  is  the  function  which  minimizes  a  quantity 
closely  related  to  the  total  curvature  of  a  curve,  subject  to  the 
constraints  that  the  spline  interpolates  the  data  and  has  two  continuous 
derivatives. 

Cubic  splines  perform  very  well  when  the  data  are  well  sampled, 
but  bathymetric  data  are  often  very  sketchy.  When  cubic  splines  were 
used  to  fit  undersampled  bathymetry  data,  the  result  was  often  something 
like  that  depicted  in  Fig.  4.  The  four  data  points  were  intended  to 
represent  the  transition  from  an  abyssal  plain  through  a  continental 
slope  to  a  continental  shelf,  as  indicated  by  the  linear  segments.  When 
the  data  are  fitted  with  a  cubic  spline,  however,  the  spline  is  seen  to 
undergo  unintended  excursions,  resulting  in  a  nonphysical  basin  and  in 
the  creation  of  an  "island".  These  anomalous  features,  hereafter  called 
spline  artifacts,  are  a  consequence  of  the  spline's  extremal i zing 
property. 


The  detection  and  correction  of  spline  artifacts  proved  to  be 
an  extremely  onerous  process,  involving  the  repeated  generation  of  high 
resolution  plots  of  suspect  portions  of  the  bathymetry  and  the 
painstaking  insertion  of  corrective  data  points  into  the  bathymetric 

data. 
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SPLINE  ARTIFACT 


FIGURE  4 

CUBIC  SPLINE  ARTIFACTS  IN  BATHYMETRY 
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The  rationale  for  the  development  to  follow  is  that,  to  be 
consonant  with  the  spline  user's  expectations,  the  T-spline  should  be 
linear  over  most  of  the  region  between  the  data  points;  nevertheless,  the 
T-spline  T(x)  is  to  be  everywhere  continuous  through  the  second 
derivative.  The  extremalizing  condition  is  counterproductive  here  and 
will  not  be  retained.  The  requirement  that  T(x)  interpolate  the  data 
will  be  slightly  relaxed,  a  concession  to  the  requirement  that  T(x)  be 
twice  differentiable  everywhere. 

To  proceed  with  the  definition  of  the  T-spline,  we  introduce 

new  abscissa  values  x.  n  and  x.  .  to  either  side  of  the  original  abscissa 

i.O  i,l 

values  xi  for  2<i<n-l  (that  is,  for  every  xi  except  the  first  one  and  the 
last  one)  as  follows: 


xi,0  =  xi  -  Ax 


(IV. 2) 


and 


x 


i.l 


x.  +  Ax 


(IV. 3) 


where  the  increment  Ax  is  always  chosen  small  enough  that 


Ax 


<  2  (xrxi-l^ 


(IV. 4a) 
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and 


(see  Fig.  4) 
new  ordinate 
lie  on  the 

<*H1,yul). 

below. 


Ax  <  2  ^Xi+l"Xi ^ 


(IV. 4b) 


.  For  each  of  these  new  abscissa  values,  we  also  introduce 
values  y.j  Q  and  y^  The  point  (xi  0,y^  q)  is  required  to 
line  segment  connecting  to  (x^  )  -  The  point 

is  required  to  lie  on  the  line  segment  connecting  (x^y^)  to 
The  T-spline  is  the  union  of  cubic  polynomials  defined 


T(x)  -  3g  +  a^ T  +  +  dgT  , 

When 

a0  =  yl-l,l* 

al  ■ 
a2  -  0, 

a3  =  0,  and 

r  =  x  -  x. _x  ^ 

T(x)  is  linear  in  this  interval;  the  cubic  polynomial 
is  degenerate. 

When  xi  g£x~xi 

a0  =  yi,0’ 

al  = 
a^  =  0 

a3  =  (yi,0"2yi+yi,i)/(fAx3),  and 
t  =  x  -  xi>Q. 
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When  x.<x<x.  , 

1  •  9 


a0  =  ^yi,0+4yi+yi,l>/6’ 
ai  =  (yi  l“yi  n)/(2Ax). 
a2  ' 

a3  =  "(^i ,o_2yi+yi ,1)/(6AX  5 »  and 

T  =  X  -  X.  . 

When  xnl<xsxi+1>0 
a0  =  yi,l’ 

al  =  (y1+l“y1 >/^xi+l“xi > • 
a2  =  0, 

a^  =  0,  and 
T  =  x'xi,l  • 

T(x)  is  linear  in  this  interval  also. 

It  is  straightforward,  although  tedious,  to  verify  that  T(x), 
thus  defined,  is  continuous  and  everywhere  twice  differentiable.  One  can 
also  show  that  T(x^)j^yi,  that  is,  T{x)  does  not  interpolate  the  data. 

C.  T-Splines  and  the  Bathymetry 

Figure  5  shows  that  when  the  bathymetry  data  of  Fig.  4  are 
fitted  with  a  T-spline,  the  "corners"  formed  by  connecting  the  bathymetry 
data  with  straight  line  segments  are  "rounded  off"  by  the  T-spline.  But, 
for  practical  bathymetries,  this  rounding  error  has  always  proven  to  be 
on  the  order  of  a  fraction  of  a  meter,  well  below  the  limits  of 

measurement  error. 

D.  T-Splines  and  Sound  Speed  Profiles 

In  general,  sound  speed  information  is  supplied  to  a  ray  model 
in  the  form  of  a  series  of  sound  speed  profiles  (sound  speeds  as 
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SURFACE  OF  OCEAN 


FIGURE  5 

T  SPLINE  FIT  TO  BATHYMETRY 
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functions  of  depth  at  given  ranges)  taken  at  various  distances  from  the 
source.  The  requirement  that  n{r,z)  and  several  of  its  gradients  be 
continuous  can  be  met  by  fitting  each  sound  speed  profile  separately  with 
T-splines,  which  provides  two  continuous  depth  derivatives  (n2  and  nzz) 
and  performing  linear  interpolations  in  range  between  profiles  to  obtain 
continuous  range  derivatives  (nr  and  nrz)  between  profiles. 

The  domains  thus  formed  are  rectangular,  except  where  the 
bottom  cuts  through  some  of  the  rectangles  to  form  its  own,  irregularly 
shaped  domain  boundaries.  Except  along  the  bottom,  the  boundaries  of  the 
rectangles  are  set  by  horizontal  lines  at  the  T-spline  node  depths  and  by 
vertical  lines  at  the  sound  speed  profile  ranges. 
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V.  EIGENRAYS 


In  principle,  it  is  possible  to  locate  an  eigenray  (a  ray  which 
connects  a  source  and  a  receiver)  by  numerically  searching  between  two 
rays  which  arrive  at  the  receiver  range  above  and  below  the  receiver 
depth.  In  practice,  it  is  not  usually  feasible  to  perform  the  numerical 
search  because  of  the  large  number  of  rays  which  would  have  to  be  traced. 
Instead,  one  can  interpolate  between  pairs  of  rays  to  locate  eigenrays. 
Also,  pairs  of  rays  suitable  for  interpolation  are  usually  not  available 
for  receivers  near  the  surface  because  it  is  unlikely  that  a  ray  will 
pass  through  the  small  depth  interval  above  the  receiver  and  below  the 
surface,  even  if  a  very  large  number  of  rays  are  traced.  On  the  other 
hand,  it  is  relatively  easy  with  a  large  sample  of  rays  to  locate  pairs 
of  rays  which  bracket  the  receiver  in  range,  so  it  is  preferable  to 
locate  eigenrays  by  interpolating  between  rays  which  bracket  the  receiver 
in  range,  not  depth. 

The  following  strategy  was  adopted  in  MEDUSA  for  finding 
eigenrays.  First,  a  large  number  of  rays  are  traced  in  small  increments 
of  the  launch  angle  from  the  source  to  a  range  beyond  the  range  of  the 
most  distant  receiver.  As  a  ray  progresses,  it  will  undergo  a  series  of 
significant  events,  such  as  surface  reflections,  transition  through 
caustics,  etc.  These  events,  their  order  of  occurrence,  and  certain  path 
information  associated  with  each  event,  together  constitute  the  ray  path 
history.  Detailed  ray  path  histories  are  recorded,  including  ranges  at 
which  rays  cross  a  receiver  depth.  After  this  ray  sweep-out  is 
completed,  eigenrays  are  located  by  scanning  the  ray  history  records  for 
pairs  of  rays  in  which  the  sequence  of  events  in  the  ray  path  histories 
are  identical,  but  which  differ  by  the  launch  angle  increment  and  which 
bracket  the  receiver  in  range  as  they  cross  the  receiver  depth.  A  third 
order  interpolation  scheme  is  used  to  determine  the  eigenray  launch 
angle.  Other  third  order  interpolations,  together  with  the  eigenray 
launch  angle  and  information  stored  in  the  history  records,  then  yield 
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arrival  angles,  travel  times,  and  path  lengths,  as  well  as  the  locations 
of  bottom  reflection  points  and  bottom  grazing  angles.  A  second  order 
interpolation  produces  eigenray  intensity.  The  details  of  the  procedure 
are  described  in  Section  V.A,  "Location  of  Eigenrays,"  and  Section  V.B, 
"Eigenray  Path  Data  Interpolation." 

This  method  of  locating  eigenrays  offers  several  advantages  in 
efficiency  and  reliability.  The  reliance  on  interpolation  reduces  the 
number  of  rays  to  be  traced  by  eliminating  the  ray  traces  which  would  be 
required  to  probe  numerically  for  eigenrays.  Comparison  of  detailed  ray 
path  histories  prevents  attempting  meaningless  interpolations  between 
rays  which  follow  significantly  different  paths  (for  example,  rays  having 
different  numbers  of  bottom  reflections).  Extrapolations  are  strictly 
avoided.  The  ray  sweep-out  approach  makes  it  possible  for  rays,  which 
may  cross  the  receiver  depth  many  times  before  the  traces  are  terminated, 
to  be  used  in  many  different  eigenray  interpolations  for  receivers  at 
different  ranges.  Ray  sweep-outs  of  several  hundred  rays  have  been  used 
to  locate  thousands  of  eigenrays. 

A.  Location  of  Eigenrays 

The  procedure  by  which  eigenrays  are  located  using  the  ray 
history  records  lends  itself  to  graphical  illustration.  Let  the  function 
r(e  )  be  the  range  at  which  a  ray,  with  launch  angle  0  ,  crosses  the 

b  b 

receiver  depth.  This  function  will  generally  be  multivalued  since  a  ray 
may  cross  the  receiver  depth  several  times  at  different  ranges.  For 
example.  Fig.  6  shows  r(9s)  for  the  sound  speed  profile  and  source  and 
receiver  shown  in  Fig.  7.  The  eigenray  launch  angles  for  a  receiver 
located  at  range  rR  are  the  zeroes  of  r(ps)-rR>  as  illustrated  in  Fig.  6 
for  rR  =  50  km. 

Devising  an  automatic  procedure  for  locating  the  zeroes 
accurately,  reliably,  and  efficiently  while  avoiding  false  predictions  is 
a  challenging  task.  As  Fig.  6  suggests,  the  range  function  r(0s)*rR  is 
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usually  a  complicated  multivalued  function  having  singularities,  slope 
discontinuities,  and  an  unknown  number  of  zeroes.  In  fact.  Fig.  6 

provides  a  comparatively  simple  illustrative  example;  a  more  typical 
range  function  will  be  discussed  in  Section  VII  (see  Fig.  13).  A 

practical  eigenray  location  scheme  must  determine  the  approximate 

locations  of  the  zeroes,  and  it  must  refine  the  launch  angle  estimates  to 
an  acceptible  degree  of  precision. 

During  the  ray  sweep-out,  r(0  )  is  sampled  at  finely  spaced 

intervals  of  o  .  If  the  sampled  points  are  designated  by  and  r., 

then,  with  some  restrictions,  the  zeroes  are  located  in  those  intervals 
where 


It  is  impractical  to  locate  an  eigenray  numerically  by  tracing 
trial  rays  in  the  interval  Osi,0s.+^1  because  of  the  time  needed  to 
compute  the  individual  rays  and  the  large  number  of  eigenrays  to  be 
found.  Instead,  r(0  )  is  replaced  in  the  interval  by  a  hermite  cubic 
interpolating  polynomial  rc(°s)>  described  below  in  Section  V.A.l,  and 
the  zero  of  r  (9  )-rD  is  located  numerically. 

As  we  shall  see,  if  hermite  splines  are  to  be  substituted  for 
r(e  ),  then  we  will  need  not  only  the  values  r.  for  each  of  the  but 

also  the  values  of  (3r/30|  )•;  that  is,  we  need  both  r  as  a  function 

Z  Zq  1 

of  0s  and  the  derivative  of  r  with  respect  to  0$  at  each  of  the  .  The 
quantity  Sr/96s I z=const  was  ^1rst  encountered  in  Chapter  II  in  connection 
with  the  calculation  of  intensity;  it  is  the  range  spreading  with 
changing  launch  angle.  We  define  £  =Dr/36s  [  z=cons^. »  much  as  we  defined 
C  =  9z/3©s  !r=const  (the  depth  spreading  function). 
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As  to  the  means  of  calculating  £,  one  might  consider  deriving  a 
differential  equation  for  £,  as  was  done  for  l,  .  That  will  not  be 
necessary,  however;  £  can  be  calculated  from  t,: 


5 


z=const 


6s=const 


r=const 


IV. 1) 


Moreover,  it  is  preferable  from  a  computational  standpoint  to  calculate 
5  from  e  rather  than  attempt  to  solve  a  differential  equation  for  £ 
because  £  is  singular  at  ray  turning  points  (where  z'=0),  as  inspection 
of  Eq.  V.l  will  quickly  disclose.  These  singularities  would  disrupt 
numerical  solution  of  the  differential  equation.  Thus,  s  serves  not  only 
in  the  calculation  of  intensity  but  also  in  the  location  of  eigenrays. 

It  turns  out  that  £'  will  also  be  needed.  Differentiation  of 
Eq.  V.l  with  respect  to  range  holding  8$  constant  yields 


C'  =  -(?'+ z"5)/z'  •  (V.2) 


1.  Hermite  Cubic  Splines 

The  hermite  cubic  spline  and  the  T-spline  have  quite 
different  applications  even  though  both  fit  data  piecewise  with  cubic 
polynomials.  Because  both  types  of  splines  are  used  extensively  and  are 
easily  confused,  the  uses  of  the  splines  will  be  distinguished  here  and 
the  hermite  cubic  spline  will  be  described. 


T-splines  are  useful  when  the  only  available  information 
about  a  function  consists  of  sample  values.  I r,  many  cases,  the  sample 
values  are  obtained  by  measurements  of  a  physical  system  and  are  expected 


to  be  representative  of  a  function  which  is  in  some  sense  smooth.  The 
T-spline  guarantees  two  continuous  derivatives  of  the  fitting  function 
and  is  used  to  fit  the  sound  speed  profiles  and  bathymetry. 

On  the  other  hand,  hermite  cubic  spline  interpolation  is 
used  to  approximate  a  known  function  when  function  evaluations  are 
costly.  It  is  more  accurate  than  linear  interpolation  because  use  is 
made  of  the  derivative  of  the  function  as  well  as  the  function  itself. 
These  splines  are  used  to  fit  the  r(e  )  curves,  for  example. 

Given  a  function  f(x),  its  derivative  f'(x),  and  the 
values  f(xQ),  f'(x0),  f(x^),  f'(x^),  the  hermite  cubic  interpolating 
polynomial  fc(x)  applicable  for  xQ£x<x^  is 


fc(x)  =  a0  +  a1(x-xQ)  +  a2(x-x0)2  +  a3(x-x0)3, 

a0=f(x0), 

ai  =  f  1  (V* 

a2  =  3[f (xx  -f(xQ)]/Ax2  -  [2f’(x0)+f1(x1)]/Ax, 
a3  =  2[f (Xq  -f(xj)]/Ax3  +  [f 1 (xQ)+f ' (x^j/Ax2,  and 


The  interpolation  error  is 


f(x)  -  fc(x)  = 


(x-xQ)  (x-Xj)2 
24  r 


(v) 


(V.3) 


for  some  v  in  [xg,x^].  This  error  is  proportional  to  Ax^  provided  f(x) 
has  four  continuous  derivatives  in  [xg,xJ.  By  contrast,  the  linear 

O  1 

interpolation  error  is  proportional  to  Ax  . 
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Derivatives  and  integrals  of  f c ( x )  can  also  be  used  to 
estimate  derivatives  and  integrals  of  f(x).  When  fc(x)  is 

differentiated,  the  interpolation  error,  f * ( x ) -f ' ( x ) ,  is  proportional  to 

1  c 
Ax  ;  differentiation  is  therefore  less  accurate  than  interpolation. 

Integration  of  f c ( x )  over  the  interval  [xg,x^]  yields  the  Euler-Maclaurin 

integration  formula  described  in  Section  III.B. 

2.  Contents  of  the  Ray  History  Records 

Each  ray  history  record  is  a  record  of  the  significant 
events  which  occurred  during  propagation  of  the  ray  in  the  order  they 
occurred,  plus  quantitative  data  about  each  event.  The  quantitative 
information  varies  with  the  nature  of  the  event.  In  this  section,  the 
events  recorded  are  described,  and  the  ancillary  data  stored  with  the 
events  is  specified.  The  reasons  for  including  specific  data  items  will 
become  apparent  when  eigenray  location  and  interpolation  procedures  are 
discussed  later  in  this  chapter.  In  general,  just  enough  information  is 
recorded  to  specify  the  location  of  the  event  and  the  ray  path  slope 
there,  to  carry  out  several  kinds  of  interpolations,  and  to  provide  ray 
path  information  which  experience  has  shown  to  be  particularly  useful. 
More  specifically,  the  ray  history  records  contain  enough  information 
about  each  ray  in  the  sweep-out  to  (1)  determine  whether  two  rays  bracket 
an  eigenray  of  a  receiver  at  a  given  range,  (2)  numerically  determine 
eigenray  launch  angles,  (3)  accurately  estimate  eigenray  arrival  angles, 
path  lengths,  travel  times,  intensities,  and  phases  and,  for  eigenrays 
which  undergo  bottom  reflections,  to  determine  the  bottom  intercepts,  ray 
path  slopes,  and  bottom  slopes  at  each  reflection  point,  and  (4)  permit 
accurate  reconstruction  of  the  ray  path  trajectories. 

a.  Surface  Reflection 


Whenever  a  ray  reflects  from  the  surface,  the  values 
of  r  and  z'  are  recorded,  along  with  z,  which  is  of  course  zero. 
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b.  Bottom  Reflection 


Upon  bottom  reflection,  the  values  of  r,  z,  z',  z", 
and  Zg  are  stored.  These  quantities  are  used  in  eigenray 
interpolations  to  determine  the  point  of  bottom  incidence  (range  and 
depth),  the  eigenray  path  slope  at  incidence,  and  the  bottom  slope  at 
incidence. 


c.  Ray  Splitting  Point 

As  shown  in  Fig.  6  and  discussed  in  Section  V.A.3, 
rays  launched  at  certain  launch  angles  may  encounter  features  in  the 
sound  speed  profile  which  cause  ray  paths  to  diverge  sharply  in  response 
to  a  slight  change  in  launch  angle,  giving  rise  to  discontinuities  in  the 
r(0s)  diagram.  The  passage  of  a  ray  through  a  point  of  divergence  is 
detectable.  If,  in  advancing  the  ray  path  solution  from  r^  to  r2,  the 
ray  path  has  an  inflection  point,  that  is,  z"  undergoes  a  sign  change, 
but  if  sign(zj-z2)  /  sign(zp,  then  the  event,  its  location,  and  the  path 
slope  there  are  recorded.  The  location  of  the  event  coincides  with  the 
inflection  point,  which  is  found  by  assuming  that  z"  varies  linearly  from 
r^  to  r2  and  finding  the  zero  of  z". 

d.  Receiver  Depth  Crossing 


i 


Whenever  a  ray  crosses  a  specific  receiver  depth,  the 
values  of  r,  z,  z',  z",  S,  S',  s,  9s/96s,  t,  and  9t/96$  are  recorded  at 
the  point  where  the  ray  crosses  the  receiver  depth.  Also  recorded  are 
the  number  of  surface  reflections,  bottom  reflections,  caustics,  turning 
points,  and  ray  splitting  points  (see  Section  V. A.2.d). 

e.  Turning  Points 

Turning  points,  where  z'=0,  are  located  by  monitoring 
the  ray  trace  for  changes  in  the  sign  of  z'.  When  one  is  detected,  the 
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approximate  range  of  the  turning  point  is  determined  by  fitting  a  hermite 
spline  to  z'  and  z"  over  the  interval  defined  by  the  last  step  of  the  ray 
trace  and  finding  the  zero  of  the  spline.  The  depth  of  the  turning  point 
is  then  obtained  by  fitting  a  hermite  cubic  spline  to  z  and  z'  and 
evaluating  the  spline  at  the  turning  point  range.  The  point  and  path 
slope  are  recorded. 

f .  Caustics 

Caustics,  points  where  ?=0,  are  located  by  monitoring 
the  ray  trace  for  changes  in  the  sign  of  C.  When  one  occurs,  the  range 
of  the  caustic  i;  located  approximately  by  fitting  a  hermite  spline  to  5 
and  V  over  the  last  range  step  and  determining  the  zero  of  the  spline. 
The  depth  of  the  caustic  is  obtained  by  fitting  a  hermite  cubic  spline  to 
z  and  z*  over  the  range  step  and  evaluating  the  cubic  for  z  at  the  range 
of  the  caustic.  The  point  and  path  slope  are  recorded. 

g.  Marker  Points 

One  of  the  most  useful  functions  of  a  ray  model  is  to 
produce  plots  of  the  ray  paths.  The  ray  history  records,  containing  as 
they  do  the  locations  of  important  events  along  the  ray  paths,  furnish 
the  information  needed  to  generate  ray  plots  and  much  else  besides.  A 
ray  path  is  plotted  by  connecting  the  events  in  its  ray  history  record 
with  a  smooth  curve  and  drawing  the  curve. 

To  take  full  advantage  of  a  ray  history  record  to 
plot  a  ray  trajectory,  one  extracts  from  the  record  the  locations  of  the 
events  in  order,  as  well  as  the  ray  path  slopes  at  the  events.  Hermite 
splines  are  used  to  connect  each  pair  of  events  in  turn,  thus 
reconstructing  the  ray  path  in  the  intervals  between  the  events.  The 
smooth  curve  formed  of  the  union  of  the  splines  is  plotted  to  complete 
the  drawing  of  the  ray  path. 


However,  there  is  no  assurance  that  significant 
events  will  occur  often  enough  to  permit  accurate  reconstruction  of  the 
ray  path.  It  is  therefore  necessary  to  insert  periodically  into  the  ray 
history  record  points  whose  only  purpose  is  to  mark  more  precisely  the 
ray  trajectory. 


In  order  to  know  when,  during  the  stepwise  solution 
of  the  ray  path  equation,  to  introduce  marker  points  into  the  history 
record,  we  need  some  means  of  determining  when  the  solution  has  advanced 
so  far  beyond  the  last  recorded  event  that  later  reconstruction  of  the 
ray  path  is  compromised.  Since  the  ray  path  is  to  be  represented  later 
by  hermite  cubic  splines,  we  need  to  know  when  the  spline  will  deviate 
unacceptably  from  the  actual  ray  path. 

Inspection  of  the  hermite  interpolation  error 
(Eq.  V.3)  suggests  that  the  largest  error  is  to  be  expected  near  the 
midpoint  of  the  interval  of  interpolation  since  the  interpolation  error 
at  the  interval  endpoints  is  zero  and  since  z1v(r)  tends  to  vary  slowly 
with  range.  But  it  is  inconvenient  and  costly  to  compute  a  new  spline 
after  each  step  and  evaluate  the  midpoint  error. 

A  much  more  satisfactory  means  of  estimating  the 
midpoint  error  is,  ironically,  to  compute  the  difference  between  the 
third  order  (cubic)  hermite  spline  and  the  much  more  complex  fifth  order 
hermite  spline  at  the  midpoint.  The  fifth  order  spline  would  use  not 
only  r,  z,  and  z',  but  also  z"  at  the  endpoints  of  an  interval.  The 
extra  information  makes  the  fifth  order  interpolation  formula  much  more 
accurate  than  the  third  order  formula,  so  that  the  difference  between  the 
third  and  fifth  order  formulas  is  a  good  estimate  of  the  error  in  the 
third  order  formula. 


But,  if  it  is  too  costly  and  inconvenient  to 
construct  the  third  order  spline  after  every  step,  then  it  is  certainly 
out  of  the  question  to  construct  both  a  third  order  spline  and  a  fifth 
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order  spline  after  every  step  in  order  to  difference  them.  However,  the 
difference  between  the  third  and  fifth  order  formulas,  when  evaluated  at 
the  midpoint,  is  a  surprisingly  simple  expression: 


*3<r0 


Ar  % 

~V  " 


z5(r 


Ar»  Ar 
2'  '  64 


2|z,(r1)-z'(r0)j  -  Ar  z"(r0)+z"(r1) 


(V.4) 


where  rg  is  the  range  of  the  last  point  entered  in  the  history  record,  r^ 
is  the  range  just  attained  with  the  last  step  of  the  ray  trace,  Ar=r^-rg, 
and  Zj  and  Zg  are  the  third  and  fifth  order  interpolation  formulas.  It 
is  not  necessary  to  construct  splines  in  order  to  make  use  of  Eq.  V.4. 
After  each  step  in  the  ray  trace,  the  right-hand  side  of  Eq.  V.4  is 
evaluated  and,  if  the  resulting  error  estimate  exceeds  a  tolerance,  the 
point  (rp*i)  is  entered  in  the  history  record  along  with  zj. 

h.  Ray  Termination 


Ray  traces  are  ended  under  several  user  controlled 
conditions.  Traces  are  terminated  if  a  ray  (1)  reaches  a  maximum 
specified  range,  (2)  exceeds  a  specified  number  of  surface  or  bottom 
reflections,  (3)  is  steeper  than  a  given  maximum  ray  path  angle,  or 
(4)  strikes  the  bottom  with  a  grazing  angle  exceeding  a  specified 
maximum.  The  ray  termination  point  and  the  path  slope  are  recorded. 

3.  Determination  of  the  Eigenray  Launch  Angle 

The  simple  r(0  )  diagram  of  Fig.  6  contains  many  of  the 
features  of  the  more  complex  diagrams  commonly  encountered  in  practice. 
Three  of  these,  described  below,  are  critical  to  the  location  of 
eigenrays  by  interpolation. 


(1)  Rays  launched  at  angles  in  the  interval  (-2.0°,  2.0°) 
never  reach  the  receiver  depth;  consequently,  no  ray 


with  a  launch  angle  in  this  interval  can  be  an 

eigenray. 

(2)  Rays  launched  at  angles  near  -5°  and  5°  encounter  a 
local  maximum  in  the  sound  speed  profile  and  diverge 
strongly.  A  discontinuity  in  r(i3$)  is  created. 

(3)  Rays  launched  at  ±10°  graze  the  surface;  rays 

launched  at  ±15°  graze  the  bottom.  At  these  angles, 
3r/3e$  is  discontinuous. 

Although  these  features  have  been  described  specifically  for  the  diagram 
of  Fig.  6,  in  general  there  are  ranges  of  values  of  0$  for  which  r(es)  is 
not  defined,  and  there  will  be  values  of  6S  for  which  r(B  )  or  3r/30s  is 
discontinuous. 


In  a  range  variable  environment  these  features  of  r(e$) 
cannot  be  so  readily  identified  with  specific  launch  angles,  yet  it  is 
not  meaningful  to  attempt  eigenray  interpolations  in  a  launch  angle 
interval  containing  any  of  these  features.  Two  rays  are  considered  to 
bracket  an  eigenray  if,  and  only  if,  (1)  they  differ  by  the  launch  angle 
increment  A0S,  (2)  they  bracket  the  receiver  in  range  when  they  cross  the 
receiver  depth,  and  (3)  the  events  in  their  ray  history  records  occur  in 
exactly  the  same  sequence,  with  the  following  exception.  If  the  history 
record  of  one  ray  shows  that  the  ray  passed  through  a  caustic  and  then  a 
turning  point,  while  the  companion  ray  passed  through  the  turning  point 
and  then  the  caustic,  it  is  nevertheless  permissible  to  use  these  two 
rays  for  eigenray  interpolation  provided  their  history  records  are 
otherwise  identical.  This  exception  is  allowed  because  no 
discontinuities  in  r(0$)  or  S(B$)  arise  from  it. 

When  two  rays  satisfy  these  requirements,  a  hermite  spline 
rc(e$)  is  computed  using  the  values  of  r^  (3r/3e$)j,  r2,  and  (3r/30s)2, 
where  the  subscripts  denote  evaluation  at  points  1  and  2  where  the  rays 
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cross  the  receiver  depth, 
by  locating  the  zero  of  rc 


The  eigenray  launch  angle  is  then  determined 
(0$)-rR  numerica1 ly. 


Successive  approximations  to  the  zero  usually  converge 
very  nearly  to  the  limits  of  precision  of  the  computer  within  three 
iterations  if  Newton's  method^8  is  used.  In  general,  if  one  is  seeking  a 
zero  of  the  function  f(x)  and  xq  is  known  to  be  a  good  approximation  to 
the  desired  zero,  then  often  an  improved  estimate,  x^,  is  given  by 
Newton's  iterative  formula 


X1  ■  x0  -  f(x0)/f‘(x0) 


(V.  5) 


To  apply  Eq.  V.5  to  the  eigenray  problem,  one  must  be  able  to  compute 
r(  9..  )-rD  and  its  derivative  dr./dO  .  But  r  (?  )  is  simply  a  cubic 

polynomial  and  drc/d6s  is  the  quadratic  polynomial  obtained  by  differen¬ 

tiating  rc(es) .  Both  functions  can  be  evaluated  very  rapidly  on  a 
computer.  Thus,  hermite  splines  and  Newton's  method  are  well  suited  to 

each  other  and  to  the  task  of  eigenray  location.  However,  it  should  be 

emphasized  that  the  true  eigenray  launch  angle  error  depends  less  on  the 
error  in  the  zero  of  the  cubic  than  the  interpolation  error  committed  by 
replacing  r(es)  with  rc(6s)»  In  order  to  keep  the  interpolation  error 
small  the  ray  sweep-out  angle  increment  A0s  must  be  kept  small. 

B.  Eigenray  Path  Data  Interpolation 

Having  determined  an  eigenray  launch  angle,  interpolations 
between  the  rays  bracketing  the  eigenray  yield  additional  information 
about  the  eigenray.  In  particular,  interpolation  gives  the  arrival 
angle,  travel  time,  path  length,  and  intensity  of  the  eigenray  when  it 
arrives  at  the  receiver.  For  eigenrays  which  undergo  bottom  reflections, 
the  reflection  points,  the  ray  path  angles  at  intercept,  and  the  bottom 
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slope  at  intercept  are  also  found.  The  interpolation  methods  are 
described  below. 

1.  Arrival  Information 


In  the  formulas  and  derivations  which  follow,  two  rays 
bracketing  the  eigenray  cross  the  receiver  depth  at  points  1  and  2,  where 
0sl<0s2*  The  ei9enray  launch  angle  is  esg. 

a.  Arrival  Angle 

A  hermite  spline  ZjUe  )  1s  defined  on  the  interval 

[0sl,0s2^  using  the  values  zi»  (3z'/3Gs)1,  z'2,  and  (az'/ae  )2-  The 
arrival  angle  is 


8r  =  tan_1(z^(0se))  .  (V.6) 


The  values  of  3z'/30s  are  not  available  directly  from  the  ray  history 
records.  They  are  computed  using  C1  as  follows: 


(V.7) 


Note  that  here  3z'/80,.  is  computed  holding  depth  constant  (z=zD),  whereas 
C  is  computed  holding  range  constant  (see  Section  II. B);  Dz'/39s  as 
computed  above  is  not  to  be  confused  with  r, '. 

The  above  expression  for  3z • /S0s  can  be  obtained  by 
resorting  to  Taylor  series  expansions.  In  Fig.  8,  two  rays  differing  in 
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launch  angle  by  5©s  (not  the  sweep-out  angle  increment  A0S)  are  shown  as 
they  pass  near  an  observation  point  (r,z).  The  ray  path  slope  at  (r,z) 
is  z1;  at  (r,z+^60s)  the  slope  is  z'  +  c'SBg  to  first  order.  The  slope  at 
(r+£Ses,z)  is  therefore  (z'+£ ' 66s)+z"60s,  again  to  first  order.  Passing 
to  the  limit  <59  -*•(), 
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(V.  8) 


where  use  has  been  made  of  Eq.  V.2. 

b.  Path  Length 

A  hermite  spline  sc  is  defined  using  the  values  s^, 
( as/aes ) J ,  s^,  and  (as/99s)2.  The  path  length  is 


s  *  sc(9 


se' 


(V.9) 
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where 


0s=const/ 


/— -j 

30 

! 

2=const/ 

=  C/cose 


(V.10) 


c.  Travel  Time 


A  hermite  spline  t  is  defined  using  the  values  t,, 
(3t/36s)p  t^,  and  (at/30  )^.  The  travel  time  is 


1  ■  ‘Ae> 

where 


n  3s 
c0  90s 


(V.ll) 


(V.  12) 


d .  Intensity 

A  hermite  spline  r  was  defined  using  r,-rD,s1,  r0-rn, 
C  I  K  |  c.  K 

and  £ 2  in  order  to  determine  the  eigenray  launch  angle: 


rc(es)  .  rc0  ♦  (es-e1)rcl  *  (9re,)2rc2  *  (e3-e,!39C3 


( V. 1 3) 
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where  rcO*  rcl>  rc2’  and  rc3  are  the  hermite  spline  coefficients.  The 
derivative  drC/,de$,  evaluated  at  0se,  is 


«c(6s>  ■  rcl  *  2(6Se-6l)rc2  +  3(ese-6l)  rc3 


( V.  14) 


Then,  by  Eq.  11.29, 
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( V. 1 5) 


2.  Bottom  Intercept  Interpolation 

In  the  following  formulas  and  derivations,  variables  with 
subscripts  1  and  2  are  evaluated  at  the  points  where  rays  intersect  the 
bottom. 


a.  Eigenray  Bottom  Intercepts 

The  range  at  which  the  eigenray  strikes  the  bottom  is 


rB  *  rBc(6se> 


(V. 16) 


where  rRc  is  the  hermite  spline  defined  by  rBl » ( 3rB/aes )  1  *  rB2’  and 
( Sr g/ 90s ) ^ •  Likewise,  the  depth  at  intercept  is 
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ZB  ~  zBc  ^  se ^ 


(V.17) 


where  zBc  is  defined  by  zgl,  (  lzg/  *  >s ) ,  zg2,  and  (fzB/o--s)2- 

An  expression  for  3rg/30  was  developed  in 

Section  II. B: 
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From  this  expression,  we  can  derive  an  expression  for  Tz.,/39  . 
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( V . 1 8) 


b.  Ray  Path  Angle  at  Intercept 

A  hermite  spline,  z|c,  is  defined  by  zjp  (3zJ/90s)^, 
z J ^ >  and  (Szj/DO^,  where  3z|/3ns  is  the  rate  of  change  of  ray  path 
slope  at  the  bottom  intercept  with  respect  to  launch  angle.  The  path 
slope  at  intercept  is 
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( V. 1 9) 
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To  find  3zj/30$  by  a  Taylor  series  expansion  we 
back  to  Fig.  3  in  Chapter  II,  where  we  see  that 


3Z| 


1  im 

&es-*o 


(V 


An  expression  for  z£,  Eg.  11.48,  was  developed  in  Section  II. B: 
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c .  Bottom  Slope  at  Intercept 

A  hermite  spline,  Zgc>  is  defined  using 
(3Zg/30s ) i ,  Zg2,  and  ( 3 Zg/S0s ) 2 *  The  bottom  slope  at  intercept  is 
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we  again  refer  to  Fig.  3  in 


To  find  DZg/30$ 


Chapter  II  and  find  that 


9ZB  _  ,,  ZB2  "  ZB1 

90s  "  68% 
s 


(V.23) 


Equation  11.49  is  an  expression  for  z' 
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ZB2  =  ZB1  +  56s  38$  ‘B1 


Combining  Eqs.  V.23  and  11.49  yields 


3zj  2j  (V.24) 
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C.  Summary 

In  this  chapter,  the  procedure  used  to  locate  eigenrays  by 
interpolation  between  rays  traced  during  a  sweep-out  was  described.  Also 
presented  were  the  interpolation  methods  used  to  determine  eigenray 
arrival  data  and  eigenray  bottom  reflection  information. 


VI.  ORGANIZATION  OF  MEDUSA 


The  ray  model  described  in  Chapters  II  through  V  of  this  reoort 
is  implemented  on  a  CDC  CYBER  171  as  a  package  of  FORTRAN  programs, 
collectively  referred  to  as  MEDUSA.  Figure  9  is  a  schematic  diagram  of 
MEDUSA,  showing  the  programs  in  the  package,  the  data  required  by  the 
programs,  and  the  output  produced  by  them. 

The  organization  of  MEDUSA  is  unusual  in  that  the  various 
functions  of  the  model  are  performed  by  separate  programs.  Each  program 
performs  a  very  limited  computational  or  display  task.  In  some  cases, 
the  output  from  one  program  is  used  as  input  by  several  other  programs. 
In  the  more  common  configuration,  a  model  is  implemented  as  a  single 
large  program,  and  all  of  the  functions  and  capabilities  of  the  model  are 
incorporated  in  the  program.  The  compartmentalized  design  of  MEDUSA, 
however,  has  several  advantages.  It  makes  the  system  more  compact  and 
reliable.  Modifications  and  additions  are  more  easily  implemented  than 
they  would  be  in  a  monolithic  program,  and  the  many  pathways  through  the 
system  make  it  possible  to  recombine  and  display  ray  information  in  a 
variety  of  ways. 

The  components  of  the  package  can  be  roughly  categorized  as 
environmental  preprocessing  and  display  routines,  ray  tracing  programs, 
and  eigenray,  propagation  loss,  and  diagnostic  postprocessors.  Their 
functions  and  relations  to  each  other  are  briefly  summarized  here.  Their 
use  is  illustrated  in  the  next  chapter. 

In  the  first  category  is  ENVPLT.  ENVPLT  plots  the  sound 
velocity  profiles  and  bathymetry.  It  is  used  for  graphical  display  of 
the  environment,  which  is  usually  received  in  the  form  of  tabulated  data. 
Such  displays  are  informative  in  their  own  right,  and  they  may  reveal 
input  data  errors. 
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Also  in  the  first  category  are  programs  which  format  data  from 
environmental  databases  for  use  by  MEDUSA.  The  author  has  assisted  in 
the  installation  of  MEDUSA  at  computer  facilities  which  possess  databases 
containing  archives  of  sound  velocity  profiles,  bathymetry,  and  bottom 
loss  data.  At  each  computer  site,  programs  convert  the  archival  data 
retrieved  by  the  database  system  into  the  format  required  by  MEDUSA. 

RAYFAN  is  the  program  responsible  for  tracing  rays  and 
generating  ray  history  records,  which  are  very  lengthy  and  so  are  stored 
on  a  disk.  In  addition  to  the  environmental  data,  RAYFAN  requires  as 
input  the  source  depth,  the  launch  angles  to  be  used  in  the  sweep-out, 
and  the  depths  of  all  receivers  for  which  eigenrays  will  be  required. 
The  ray  history  records  produced  are  used  by  a  growing  array  of 
postprocessors. 

RAYTREK  plots  the  ray  paths  of  any  of  the  rays  recorded  in  the 
ray  history  records. 

CAUSTIC  plots  the  locations  of  caustics  encountered  during  the 
ray  sweep-out. 

THSPLT  produces  plots  of  r(e$);  THRPLT  plots  arrival  angles  as 
a  function  of  receiver  range. 

SERPENT  locates  eigenrays  by  interpolation  of  the  ray  histories 
and  generates  eigenray  path  history  records.  The  eigenray  records  are 
also  extremely  lengthy  and  therefore  are  stored  on  a  disk  also.  EIGONY 
prints  eigenray  path  summaries  for  specified  receivers. 

PROPLOS  adds  the  pressure  contributions  of  the  eigenrays  and 
computes  propagation  loss.  PLPLOT  plots  the  propagation  loss  curves. 
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This  chapter  illustrates  the  use  of  the  MEDUSA  programs.  For 
the  example  presented  here,  the  environment,  particularly  the  bathymetry, 
changes  rapidly  with  range,  bottom  interaction  effects  are  important,  and 
the  track  ranges  are  relatively  short.  MEDUSA  has  also  been  used  in  long 
range  propagation  analysis  problems.  The  computer  resources  required 
have  proved  to  be  remarkably  insensitive  to  the  ranges  involved.  For  the 
long  range  case,  long  ray  traces  are  compensated  by  the  reduced  number  of 
rays  to  be  traced  and  the  small  path  angles  involved. 

Each  use  of  MEDUSA  begins  with  the  generation  of  environmental 
plots  by  ENVPLT.  Figure  10  is  the  graphical  output  from  ENVPLT  for  this 
example.  The  bathymetry  and  sound  speed  profiles  are  presented  to  the 
programs  as  tabulated  data  in  the  usual  fashion.  The  tabulated  sound 
speed  and  bathymetry  points  supplied  to  ENVPLT  are  plotted  as  discrete 
points.  The  smooth  curves  interpolating  the  user  supplied  data  points 
are  generated  by  repeated  evaluation  of  the  bathymetric  and  sound  speed 
profile  T-splines  at  finely  spaced  intervals  of  depth  and  range.  Several 
sound  speed  profiles  are  also  shown  in  Fig.  10  between  the  user  supplied 
profiles.  These  profiles,  distinguished  by  the  absence  of  accompanying 
discrete  data  points,  are  intermediate  profiles  which  show  how  the 
profiles  evolve  with  range  from  one  user  supplied  profile  to  the  next. 
Environmental  plots  thus  provide  not  only  a  graphical  representation  of 
the  somewhat  cryptic  tabulated  data,  but  are  also  a  valuable  diagnostic 
aid  against  input  data  errors,  undersampling  of  the  environment,  and 
possible  range  interpolation  anomalies  in  the  sound  speed  profiles. 

Occasionally,  usually  for  special  research  purposes,  a  user 
will  supply  sound  speeds  or  bathymetry  to  MEDUSA  as  mathematical 
functions  of  his  own  devising,  bypassing  the  standard  procedure  of 
fitting  T-splines  to  tabulated  data.  ENVPLT  will  then  produce  plots  of 
the  user  defined  environmental  functions.  The  substitutions  pose  no 
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FIGURE  10 

SOUND  SPEED  PROFILES  AND  BATHYMETRY  (ENVPLT) 
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special  problems  for  MEDUSA,  provided  the  functions  possess  the  required 
degrees  of  continuity  discussed  in  Chapters  III  and  IV. 

Once  satisfied  with  MEDUSA'S  representation  of  the  environment, 
the  next  step  is  usually  to  produce  ray  path  plots,  such  as  Fig.  11.  In 
this  example,  the  source  is  located  1  m  above  the  bottom  at  a  depth  of 
5303  m.  This  is  accomplished  by  tracing  a  few  rays  with  RAYFAN  and 
plotting  the  ray  paths  using  RAYTREK.  These  plots  are  often  all  that  is 
required  of  MEDUSA,  but  they  can  also  convey  very  quickly  to  an 
experienced  user  a  great  deal  of  information  about  the  effects  of 
bathymetry  and  range  and  depth  variations  of  the  sound  speed  on  sound 
propagation.  Thus,  RAYTREK  is  often  used  as  a  survey  tool  in  preparation 
for  eigenray  and  propagation  loss  calculations. 

In  order  to  proceed  beyond  the  generation  of  ray  path  plots,  a 
ray  sweep-out  is  performed  with  RAYFAN,  in  which  hundreds  of  rays  are 
traced  in  fine  increments  of  the  launch  angle  to  a  maximum  range  beyond 
the  range  of  the  most  distant  receiver  of  interest.  The  ray  history 
records  generated  form  the  basis  for  all  further  processing.  In  this 
example,  the  rays  traced  had  launch  angles  ranging  from  -50°  to  50°  in 
0.5°  increments,  except  for  the  interval  from  -10°  to  10°,  where  a  finer 
0.25°  increment  was  used.  A  receiver  depth  of  18.3  m  (60  ft)  was 

specified.  Rays  were  traced  to  a  maximum  range  of  35  km. 

Figure  12,  generated  by  program  CAUSTIC,  shows  the  caustics 

encountered  by  the  rays  during  the  sweep-out.  Unusually  high  signal 
pressure  levels  occur  near  caustics,  which  are  focal  points  of  acoustic 
energy.  Although  the  ray  theory  approximation  K'k  is  invalid  within  a 

few  wavelengths  of  a  caustic,  the  effect  of  the  approximation  is  always 
to  cause  uncorrected  ray  theory  to  overestimate  the  pressure  near  a 
caustic.  Thus,  ray  theory  predictions  of  the  locations  of  acoustic  focal 
points  are  accurate,  although  predictions  of  actual  signal  pressure  are 

exaggerated  extremely  close  to  them.  Corrections  can  be  made  to  revise 
ray  theory  pressure  predictions  downward  and  such  corrections  are  likely 
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to  be  introduced  into  MFDUSA  at  a  later  time.  In  any  case,  Fig.  12  shows 
the  caustics  are  well  below  the  18.3  m  receiver  depth  and  we  need  not  be 
concerned  with  them  in  this  example. 


THSPLT  produced  the  r(  )  plot  of  Fig.  13.  Each  point  on  the 
plot  represents  a  range  at  which  a  ray  crossed  the  receiver  depth  and  the 
launch  angle  of  the  ray.  The  fine  sampling  of  launch  angles  causes  the 
points  to  pack  closely  together,  so  that  they  tend  to  form  smooth  curves, 
although  the  curves  do  have  discontinuities  and  cusps.  The  intersections 
with  the  curves  of  a  straight  line,  drawn  across  the  plot  at  a  constant 
range  corresponding  to  a  receiver  range,  occur  at  launch  angles 
corresponding  to  eigenray  launch  angles.  By  moving  the  straight  line 
back  and  forth  over  the  plot  to  represent  receivers  at  various  ranges, 
one  can  observe  the  eigenray  launch  angles  changing  as  the  line 
intersects  different  parts  of  the  curves. 

The  arrival  angle  diagram  of  Fig.  14  was  produced  by  THRPLT  in 
much  the  same  way  as  THSPLT  generated  the  launch  angle  diagram  of 
Fig.  13;  the  points  on  the  arrival  angle  diagram  indicate  the  ranges 
where  rays  crossed  the  receiver  depth  and  the  angles  the  ray  paths  made 
with  respect  to  the  horizontal  as  they  crossed  the  receiver  depth.  The 
intersections  of  the  arrival  angle  curves  with  a  straight  line,  drawn 
across  the  plot  at  a  given  receiver  range,  indicate  the  expected  ray 
arrival  angles  for  a  receiver  at  that  range.  By  varying  the  receiver 
range,  one  can  tell  at  a  glance  which  arrival  angles  to  expect  at  various 
receiver  ranges  and  how  the  arrival  angle  structure  may  be  expected  to 
change  with  range. 

The  preceding  display  programs,  except  ENVPLT,  all  simply 
plotted  information  stored  in  the  ray  history  files.  To  proceed  beyond 
this  point,  further  processing  is  necessary.  The  next  step  is  to  locate 
eigenrays,  using  SERPENT  to  perform  eigenray  interpolations  on  the  data 
stored  in  the  ray  history  records.  SERPENT  will  produce  eigenray  history 
records. 
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r(0s)  FOR  RANGE  DEPENDENT  MEDIUM  (THSPLT) 
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ANGLES  AT  MOVING  SENSOR  -  deg 


EIGONY  will  display  any  desired  portion  of  the  eigenray 
records.  Figure  15  shows  a  very  brief  excerpt  from  the  output  normally 
generated  by  EIGONY.  Here,  we  see  that  SERPENT  found  eight  eigenrays  for 
the  receiver  located  at  a  range  of  14  km  and  a  depth  of  18.3  m.  For  each 
of  those  eigenrays,  EIGONY  prints  the  launch  and  arrival  angle,  the 
travel  time,  the  number  of  surface  and  bottom  reflections  the  eigenray 
underwent  and  the  number  of  caustics  it  passed  through  on  the  way  to  the 
receiver,  and  the  acoustic  pressure  at  the  receiver.  The  pressure  is 
computed  according  to  the  geometrical  spreading  law  of  Eq.  11.29;  the 
frequency  dependent  effects  of  volumetric  attenuation  in  the  water  and 
bottom  loss  are  accounted  for  later.  For  each  bottom  reflection  which  an 
eigenray  undergoes,  EIGONY  prints  the  location  of  the  bottom  intercept, 
the  bottom  slope  angle,  and  the  grazing  angle  (the  angle  of  the  ray  with 
respect  to  the  bottom  at  the  intercept). 

PROPLOS  adds  the  complex  eigenray  pressure  contributions  at 
each  of  the  receiver  locations.  If  the  user  specifies  tables  of 
reflection  coefficients  (which  are  generally  frequency  dependent),  then 
the  reflection  loss  and  phase  shift  due  to  each  bottom  reflection  of  each 
eigenray  are  figured  into  the  pressure  summation.  At  the  same  time,  if 
the  user  wishes  to  apply  a  frequency  dependent  volumetric  attenuation 
coefficient,  then  those  losses  are  taken  into  account  also.  The  total 
complex  pressure  P  at  a  given  receiver  is 


P  = 
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(VII. 1) 


where  N  is  the  number  of  eigenrays,  a  is  the  attenuation  (dB/m),  sm  is 
the  path  length  of  the  mth  eigenray,  tm  is  the  travel  time  of  the  mth 
eigenray,  n  is  the  number  of  surface  reflections,  n  is  the  number  of 
caustics,  Nb  is  the  number  of  bottom  reflections,  and  R^  is  the  complex 
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reflection  coefficient  for  the  Uh  bottom  reflection  of  the  mth  eigenray. 
The  phase  inversion  of  a  ray  upon  surface  reflection  arises  from  the 
representation  of  the  surface  as  a  flat  pressure  release  perfect 
reflector.  Notice  the  -tt/2  phase  shift  as  a  ray  passes  through  a 
caustic.8  The  sign  of  the  phase  shift  is  often  reversed  (incorrectly)  in 
ray  models,  and  the  phase  shift  itself  is  sometimes  mistakenly  associated 
with  turning  points  rather  than  caustics. 

The  coherent  propagation  loss  is 
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The  incoherent  propagation  loss  is  defined  to  be 
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Here,  the  root-mean-square  pressures  are  added  without  regard  to  phase. 
This  procedure  is  not  well  justified  mathematically,  but  it  does  yield 
propagation  loss  curves  (the  "incoherent"  propagation  loss)  which  vary 
more  slowly  with  range  than  the  coherent  propagation  loss,  and  which,  in 
some  cases,  seem  to  approximate  a  range  averaged  coherent  propagation 
loss.  In  any  event,  the  incoherent  propagation  loss  is  frequently 
required  of  a  ray  model  and  must  be  provided. 

PROPLOS  computes,  prints,  and  stores  on  disk  the  coherent 
propagation  loss,  the  incoherent  propagation  loss,  and  the  incoherent 
propagation  loss  computed  assuming  an  infinitely  rigid,  perfectly 
reflecting  bottom. 
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PLPLOT  plots  the  propagation  loss  computed  by  PROPLOS. 
Figure  16  shows  the  coherent  and  incoherent  propagation  loss  curves  for  a 
frequency  of  150  Hz  and  an  attenuation  of  1.5  x  10"®  dB/m.  Bottom  loss 
tables  were  also  specified,  so  that  each  time  an  eigenray  struck  the 
bottom  it  incurred  some  grazing  angle  dependent  loss  but  no  phase  shift. 
Figure  16  also  shows  the  incoherent  propagation  loss  for  an  infinitely 
rigid,  perfectly  reflecting  bottom.  The  difference  between  the 
incoherent  propagation  loss  curves  for  the  reflecting  bottom  and  the 
partially  absorbing  bottom  gives  a  quick  indication  of  the  importance  of 
the  bottom  in  propagation  analysis. 
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SOURCE  DEPTH  =  5303.000  m,  RECEIVER  DEPTH  =  18.300  m 
FREQUENCY  =  150.000  Hz,  ATTENUATION  =  0.000010  dB/m/kHz 


COHERENT  PROPAGATION  LOSS 
INCOHERENT  PROPAGATION  LOSS 


- INCOHERENT  (RIGID  BOTTOM) 


FIGURE  16 

PROPAGATION  LOSS  (PLPLOT) 
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VIII.  CONCLUSION 


When  ray  models  began  to  make  the  transition  from  range 
invariant  to  range  dependent  media,  many  of  the  techniques  developed  for 
range  independent  media  were  no  longer  applicable.  In  range  independent 
environments,  the  ray  paths  are  periodic  functions  of  range,  and  a  ray 
with  a  given  launch  angle  has  a  characteristic  period,  or  cycle  distance. 
The  first  ray  models  fully  exploited  that  fact  to  greatly  simplify  ray 
tracing  and  locating  eigenrays  (see  Ref.  11). 

Since  rays  in  a  horizontally  stratified  medium  are  cyclic,  it 
is  only  necessary  to  perform  the  calculations  needed  to  trace  a  ray 
through  half  of  a  cycle;  the  ray  trajectory  can  then  be  extrapolated  to 
any  desired  range  by  connecting  half  cycle  paths  and  segments  of  half 
cycle  paths.  Thus,  range  independent  ray  models  could  trace  rays  very 
rapidly.  Most  such  models  used  analytical  solutions  of  the  ray  equations 
(see  Section  III. A),  although  this  was  not  necessary  and,  in  retrospect, 
probably  set  an  unfortunate  precedent. 

It  was  possible  and  very  useful  in  horizontally  stratified 
media  to  be  able  to  classify  rays  and  eigenrays  according  to  just  four 
criteria:  the  sign  of  the  ray  launch  angle  (whether  the  ray  was  upgoing 
or  downgoing  at  the  source),  the  sign  of  the  arrival  angle  (whether  the 
ray  was  upgoing  or  downgoing  at  the  receiver),  the  number  of  cycles,  and 
whether  the  ray  reflected  from  the  surface  or  bottom  or  had  turning 
points  instead  of  relections.  The  launch  angle  intervals  giving  rise  to 
the  various  ray  categories  could  be  readily  determined  from  the  sound 
speed  profile,  Snell's  law,  and  the  range.  It  was  also  possible  to 
determine  whether  there  were  any  eigenrays  belonging  to  a  particular 
classification  for  a  receiver  at  a  particular  range.  Thus,  eigenrays 
could  be  located  by  systematic  exploration  of  each  of  the  ray  categories. 
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It  was  even  possible  to  locate  an  eigenray  of  known  classification  by 
numerical  shooting  methods  (see  Ref.  11). 

But  rays  in  range  variable  environments  are  not  periodic 
functions,  and  Snell's  law  is  no  longer  applicable.  New  strategies  for 
tracing  rays  and  locating  eigenrays  in  these  more  challenging 
environments  have  been  the  subject  of  this  report. 

It  is  no  longer  sufficient  in  such  media  to  perform  ray  trace 
calculations  over  a  half  cycle;  rays  must  be  traced  out  in  full,  at 
substantially  greater  computational  cost.  There  are  compelling  reasons 
for  seeking  numerical,  rather  than  analytical,  solutions  to  the  ray 
equations.  With  that  in  mind,  the  equations  for  the  ray  path  trajectory 
and  intensity  are  written  in  terms  cf  differential  equations  in  a  form 

1  Q 

proposed  by  Solomon  and  Armijo,  because  a  great  variety  of  methods 
exist  for  solving  differential  equations  numerically.  Chapter  III 

discusses  the  numerical  alternatives  in  terms  of  the  particular  problem 
at  hand  and  arrives  at  a  suitable  method  which,  although  somewhat  slow, 
is  very  accurate,  flexible,  simple,  and  reliable.  The  method,  Kutta's 

third  order  rule,  lacks  a  procedure  for  estimating  errors.  An 
independent  means  of  error  estimation  and  control  is  therefore  developed. 

Some  of  the  ray  path  information,  namely  the  path  length  and 

travel  time,  can  be  conveniently  expressed  as  ray  path  integrals  and 

computed  numerically.  An  appropriate  numerical  integration  method,  based 
on  an  Euler-Maclaurin  formula,  is  explained. 

When  rays  reflect  from  the  surface  or  bottom,  the  solutions  of 
the  ray  path  and  intensity  differential  equations  must  be  restarted  at 

the  boundary  intersection  with  new  initial  values.  The  reinitialization 
conditions  for  the  ray  path  equation  are  easily  derived  from  geometrical 
considerations  since  it  is  assumed  that  rays  undergo  specular  reflection 
from  boundaries.  But  the  effects  of  reflections  from  boundaries  on  ray 
path  spreading  are  more  complicated.  The  reinitialization  equations  for 
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ray  path  spreading,  derived  in  Section  II.B,  correctly  account  for  the 
effects  of  ray  path  curvature,  and  for  the  curvature  of  the  bottom  when 
applicable. 

The  solution  of  the  ray  equations  is  very  sensitive  to  the 
representation  of  the  sound  speed  and  bathymetry  within  the  model,  and 
the  location  of  eigenrays  depends  critically  on  a  condign  treatment  of 

the  environment.  Sound  speeds  and  bathymetry  are  normally  supplied  to  a 

ray  model  as  tables  of  values  measured  at  various  ranges  and  depths.  The 

ray  model  must  fit  the  tabulated  data  with  interpolating  functions  so 

that  the  sound  speed  and  bottom  depth  can  be  determined  away  from  the 
tabulated  points.  The  interpolating  functions  should  possess  two 
continuous  derivatives  in  order  to  avoid  propagation  loss  anomalies  and 
to  ease  the  task  of  eigenray  location.  The  functions  and  their 

derivatives  should  be  easy  to  evaluate  on  a  computer. 

Cubic  splines  had  been  proposed  as  suitable  interpolating 

functions,  and  had  even  been  incorporated  into  some  ray  models.  But  when 
the  author  applied  cubic  splines  to  tabulated  bathymetry,  it  was  soon 

discovered  that  the  splines  would  usually  exhibit  nonphysical 

contortions.  Attempts  to  suppress  the  formation  of  these  spline 

artifacts,  by  tediously  inserting  new  points  into  the  bathymetry  tables, 
often  ended  in  frustration  when  the  spline  would  perversely  introduce  new 
artifacts  elsewhere.  Spline  artifacts  also  appeared  occasionally  in 
sound  speed  profiles. 

The  problem  of  spline  artifacts  was  eventually  solved  by 
abolishing  the  cubic  splines  which  gave  rise  to  them  in  favor  of  modified 
splines  devised  by  the  author,  called  T-splines.  T-splines  and  the 
environment  are  discussed  in  Chapter  IV. 

One  of  the  principal  duties  of  a  ray  model  is  to  locate 
eigenrays  and  determine  the  ray  path  histories  of  the  eigenrays.  With 
this  capability,  the  user  of  a  ray  model  can  identify  the  paths  along 
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which  energy  is  delivered  to  the  receiver.  Eigenrays  are  located  in  a 
range  dependent  medium  by  interpolating  between  rays*  traced  during  an 
extensive  ray  sweep-out,  which  bracket  a  receiver  in  some  sense.  The 
precise  determination  of  when  it  is  appropriate  to  make  eigenray 
interpolations  and  how  to  obtain  eigenray  path  information  by 
interpolation  are  two  of  the  major  achievements  discussed  in  this  thesis 
(see  Chapter  VI).  Ray  history  records  of  sequences  of  hundreds  of 
judiciously  selected  events,  such  as  bottom  reflections  and  passage 
through  caustics,  contain  information  crucial  to  both  enterprises.  The 
concept  of  the  ray  history  record  derives,  greatly  elaborated,  from  the 
ray  classification  scheme  of  range  independent  models,  which  had  only 
four  elements.  Eigenray  path  data  interpolation  methods  make  use  of 
hermite  cubic  interpolator  splines  to  fit  data  stored  in  the  history 
records,  and  are  far  more  accurate  than  the  first  order  finite  difference 
methods  commonly  used. 


Ray  models  traditionally  have  played  an  important  role  in 
underwater  acoustics  research,  and  it  is  anticipated  that  their  role  will 
expand  as  their  capabilities  increase  and  new  applications  are  found  for 
them.  Some  of  the  improvements  and  extensions  likely  to  be  undertaken 
soon  are  described  below. 


It  is  expected  that  efforts  to  develop  significantly  faster  ray 
trace  methods  will  be  made  in  the  near  future.  Success  in  this  rather 
prosaic  matter  would  nevertheless  greatly  encourage  the  distribution  and 
usage  of  ray  models. 

Range  dependent  ray  models  are  poised  to  receive  a  plethora  of 
wave  theory  modifications  and  corrections  which  are  designed  to  improve 
their  treatment  of  caustics,  ducts,  shadow  zones,  and  surface  and  bottom 
interactions.  Some  range  invariant  ray  models  already  contain  these 
improvements;  it  remains  only  to  generalize  them  to  range  dependent 
media. 
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The  MEDUSA  model  is,  at  present,  restricted  to  cylindrically 
symmetrical  environments  but,  with  minor  modifications,  it  will  be 
extended  to  certain  kinds  of  three-dimensional  environments  which  possess 
some  form  of  azimuthal  or  translational  symmetry.  The  preliminary 
studies  of  propagation  across  idealized  eddies,  fronts,  and  seamounts, 
and  of  oblique  propagation  over  ridges  and  continental  slopes  can  be 
conducted  without  resort  to  fully  three-dimensional  models.  For  this 
reason,  the  development  of  general  purpose  three-dimensional  models  may 
be  premature  at  present. 

It  is  hoped  that  this  report  will  serve  to  describe  some  of  the 
practical  matters  which  must  be  considered  when  implementing  or  simply 
using  ray  models  (particularly  the  MEDUSA  ray  model),  and  to  suggest 
areas  of  further  research  which  would  extend  the  capabilities  and 
applications  of  ray  models. 
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